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SECTION  I 


INTRODUCTION 


1.1  THE  PROBLEM 

Use  of  filamentary  composites  in  aircraft  components  has  been  steadily  growing. 
For  satisfactory  design  it  is  imperative  that  the  service  behavior  and  the  ultimate  load 
capacity  of  the  components  be  well  understood.  Experimental  procedures  cannot 
realistically  be  expected  to  cover  all  possible  stress  configurations  for  all  possible 
designs.  Analytical  models  of  material  behavior,  which  can  be  used  to  predict  structural 
behavior  under  a  variety  of  service  conditions,  are  necessary  to  develop  a  firm  basis 
for  design. 

1.2  RESEARCH  OUTLINE 

The  purpose  of  the  research  described  in  this  report  was  to  develop  reliable  and 
efficient  computational  procedures  for  prediction  of  deformation  and  stress  causing 
cumulative  damage  in  structural  components  made  of  composite  laminates.  The  proposed 
models  were  to  be  verified  in  an  appropriate  manner. 

The  w’ork  performed  under  the  research  program  consisted  primarily  of  analytical 
studies  involving  mathematical  modelling  of  behavior  of  laminated  plates  under 
combined  bending  and  stretching  and  of  free-edge  delamination  specimens  subjected  to 
uniform  extension. 

The  models  developed  were  designed  to  be  capable  of  representing  simultaneous 
occurrence  of  flexure  and  stretching  and  have  the  capability  of  predicting  interlaminar 
norma!  and  shear  stresses. 


1 


Review  of  available  literature  pointed  out  the  inadequacy  of  the  existing  theories 
of  mechanics  of  laminated  composites.  Theories  based  on  assumption  of  a  linear  or 
higher  order  variation  of  the  in-plane  and  the  transverse  components  of  displacement 
over  the  thickness  of  the  laminate  have  been  known  to  be  incapable  of  representing 
local  effects.  Even  the  ’’layerwise”  or  discrete  laminate  theory  could  not  properly 
predict  interlayer  stresses  and  there  was  no  satisfactory  procedure  to  allow  for  shear 
deformation  effects.  Pagano's  [1978]  assumed  stress  distribution  approach  is  theoretically 
sound  but  its  use  to  find  numerical  solutions  to  specific  problems  was  limited  to  only 
a  few  layers.  The  finite  element  methods  available  could  not  accurately  model  the 
complete  stress  state  especially  at  interfaces  and  free  edges. 

The  research  reported  herein  attempts  to  plug  some  of  the  gaps.  An  hierarchical 

approach  to  the  construction  of  a  sequence  of  refined  discrete  laminate  theories  to 

model  combined  bending  and  stretching  of  laminated  plates,  based  upon  assumed 
displacement  patterns  was  developed.  Coupled  constitutive  equations  for  force  resultants 
were  derived  using  an  extension  of  Reissner's  variational  procedure.  This  eliminates  the 
need  for  ad  hoc  ’’correction  factors”  to  allow  properly  for  shear.  Pagano's  theory  was 
restated  in  a  self-adjoint  form  with  a  reduction  in  the  number  of  free  field  variables. 
This  renders  the  approach  convenient  for  variational  formulation  and  application  of 
finite  element  techniques.  For  analysis  of  free-edge  delamination  specimens,  continuous 
traction  finite  elements  were  developed.  Test  programs,  for  which  adequate 

documentation  was  available,  were  carefully  examined.  The  predicted  modes  of  failure 
and  progressive  damage  sequence  were  compared  with  the  test  data.  Figure  1  shows 
the  various  components  of  the  research  program  and  the  individual  reports,  listed  in 
Appendix  B,  wherein  each  is  fully  documented. 


2 


C/1 


3 


Fig  ire  1:  SUMMARY  OF  THE  RESEARCH  EFFORT. 


1.3  STRUCTURE  OF  THE  REPORT 

This  report  describes  the  overall  research  effort  and  the  principal  results  of  the 

investigation.  We  first  present,  in  Section  II,  a  summary  of  previous  work  on 

modelling  behavior  of  composite  laminates.  Section  III  contains  a  summary  of  the  basic 
equations  of  linear  elasticity  and  their  specialization  to  the  case  of  free-edge 

delamination  specimens.  Section  IV  presents  an  hierarchical  approach  to  theory  of 

laminated  plates  using  assumed  displacement  patterns.  Section  V  gives  a  summary  of 
the  new  discrete  laminate  theories,  developed  in  the  present  research  program,  allowing 
for  coupled  constitutive  equations  for  the  force  resultants.  Section  VI  contains  a 
restatement  of  Pagano’s  assumed  stress  distribution  theory.  Section  VII  points  out  the 

self-adjoint  form,  useful  for  variational  formulation,  of  laminate  theories.  The  governing 

equations,  for  the  collection  of  laminate  theories,  including  the  field  and  interface 

continuity  equations  are  stated  in  a  self-adjoint  form  and  general  variational  principles 
are  proposed.  Extensions  of  the  general  formulation,  to  admit  relaxation  of  continuity 
requirements  on  some  of  the  field  variables,  are  indicated  as  well  as  specializations  to 
reduce  the  number  of  free  field  variables  by  requiring  some  of  them  to  identically 
satisfy  some  of  the  field  equations.  This  is  often  necessary  to  make  the  problem 
computationally  tractable.  Section  VIII  presents  a  summary  of  finite  element  studies. 
Appendix  A  contains  a  summary  of  the  general  procedure  for  setting  up  variational 

formulations  for  linear  coupled  problems.  Appendix  B  lists  the  publications  that  have 
resulted  from  the  research  effort. 


SECTION  II 


REVIEW  OF  PREVIOUS  WORK 

2.1  INTRODUCTION 

This  section  briefly  reviews  previous  work  on  modelling  behavior  of  composite 
laminates  with  a  view  to  identifying  cumulative  damage  processes.  Two-dimensional 
representation  of  the  three-dimensional  physical  problem,  and  finite  element  solution 
procedures  applied  to  evaluation  of  stress  and  deformation  in  composite  laminates  are 
discussed.  Detailed  reviews  are  available  in  the  literature.  As  examples  one  may  cite 
the  recent  reviews  by  Reddy  [1990]  and  Hanna  [ 1 990].  The  former  is  in  the  nature  of 
a  survey  while  the  latter  presents  details  of  various  theories  along  with  applications. 

22  TWO-DIMENSIONAL  REPRESENTATION 

2.2.1  Introduction 

In  order  to  simplify  solution  of  the  three-dimensional  problem  of  stresses  and 
deformations  in  composite  laminates,  it  is  desirable  to  take  note  of  certain  material  and 
geometric  symmetries  that  might  exist  and  the  fact  that  in  most  applications  the 
thickness  of  the  laminate  is  very  small  in  comparison  with  the  lateral  dimensions. 
For  homogeneous  plates,  it  is  customary  to  make  the  dependence  of  the  field  variables 
upon  the  transverse  coordinate  explicit.  An  assumption  is  usually  made  for  the  pattern 
of  variation  of  the  displacement  or  the  stress  components  with  this  coordinate.  For  the 
case  of  free-edge  delamination  specimens  under  uniform  axial  strain,  the 
three-dimensional  problem  reduces  to  a  pseudo-two-dimensional  one  because  of  the 
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dependence  upon  the  longitudinal  coordinate  x,  being  eliminated.  Specialization  of  the 
thin  plate  theories  to  the  case  of  free-edge  delamination  specimens  makes  the  problem 
one-dimensional. 

2.2.2  Displacement-based  Theories. 

Existing  displacement-based  theories,  as  applied  to  laminated  composite  plates,  can  be 
classified  into  three  categories,  viz., 

1.  Classical  thin  plate  theory  and  the  first  order  shear  deformation  theory. 

2.  Higher  order  theories. 

3.  Discrete  laminate  theories. 

The  first  two  categories  assume  the  dependence  of  the  components  of  the 
displacement  vector  to  be  polynomials  in  the  coordinate  normal  to  the  surface  of  the 
plate.  The  third  category  treats  each  lamina  as  a  homogeneous  anisotropic  plate  and 
enforces  continuity  of  components  of  traction  and  displacement  at  the  interfaces. 

2.2.2. 1  Classical  Thin  Plate  Theory  [CPTl 

A  theory  for  a  laminate  was  presented  by  Reissner  and  Stavsky  [1961]  for  a 
two-layer  system.  Stavsky  [1961]  extended  it  to  a  multi-layer  plate.  Dong  et  al.  [1962] 
extended  the  concepts  to  the  analysis  of  anisotropic  laminated  shells.  This  theory 
assumes  the  in-plane  components  of  displacement  to  be  linear  in  the  transverse 
coordinate  and  the  out-of-plane  component  to  be  independent  of  it.  Also,  the  rotation  of 
any  cross-section  is  assumed  to  be  equal  to  the  rotation  of  the  tangent  at  the 
midsurface.  This  is  equivalent  to  ignoring  shear  deformation  and  assuming  the  thickness 
of  the  plate  to  remain  constant.  The  transverse  strain  is  zero  by  the  displacement 
asumptions  and  the  t.  ansverse  stress  is  also  assumed  zero.  Implying  infinite  shear 
rigidity,  it  leads  to  an  overestimation  of  plate  stiffness  resulting  in  underestimation  of 
transverse  deflection  [Whitney,  1969;  Pagano,  1969,  1970;  Srinivas  and  Rao,  1970].  The 
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effect  is  more  pronounced  in  composite  laminates  due  to  the  low  ratio  of  shear 

modulus  to  Young's  modulus.  Also,  the  theory  cannot  predict  local  and  interfacial 

behavior  of  the  laminate  which  is  critical  to  effectively  model  damage  associated  with 

delamination. 

2.2.22  First  Order  Shear  Deformation  Theory  [FSDT]. 

In  the  first  order  shear  deformation  theory  [e.g.,  Yang  et  al.,  1966;  Whitney  and 
Pagano,  1970]  the  rotation  of  the  cross-sections  is  allowed  to  be  different  from  the 
rotation  of  the  tangent  at  the  midsurface.  This  is  an  extension  of  Reissner's  [  1 944, 
1945]  and  Mindlin's  [  1 951  ]  ideas.  It  can  allow  for  shear  deformation  and  predict  lateral 
deflections  and  fundamental  frequencies  of  vibration  reliably  but  does  not  give  good 
quantitative  results  [Whitney  1969;  Srinivas  and  Rao,  1970]  for  the  transverse  shear 
stress  distribution.  According  to  Lo  et  al.  [1977],  ’’Despite  the  increased  generality  of  the 
shear  deformation  theory,  the  related  flexural  stress  distributions  show  little 
improvement  over  these  of  the  classical  plate  theory.”  To  improve  accuracy,  Whitney 
and  Pagano  [1970]  and  Kulkarni  and  Pagano  [1972]  suggested  use  of  a  suitable  shear 
correction  factor  K.  However,  the  value  of  this  factor  was  not  a  fixed  constant 
[Whitney  and  Pagano  1970;  Kulkarni  and  Pagano  1972]  and  different  values  w’ere  seen 
to  be  optimal  for  different  layups  and  for  various  quantities  of  interest  in  the 

particular  problem.  Whitney  and  Pagano  [1970]  stated;  ’The  evaluation  of  K  in  a 

specific  problem  depends  on  either  the  exact  elasticity  solution  of  the  problem  or 
experimental  evidence.  Owing  to  the  scarcity  of  such  information,  we  can  say  little 
about  the  evaluation  of  K  at  the  present  time.”  In  Mindlin's  approach  [1951],  thickness 
shear  modes  obtained  by  the  approximate  theory  and  by  an  exact  elasticity  solution 

for  straight-crested  flexural  waves  in  an  infinite  plate  w-ere  matched.  Yang  et  al. 
[1966]  applied  the  general  theory  to  propagation  of  plane  strain  waves  of 

Rayleigh-I^imb  type  in  specific  two-layer  isotropic  plate  and  showed  that  the  value  of 
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k  depends  upon  the  material  properties  of  the  layered  plate.  Nelson  and  Lorch  [1974] 
presented  a  refined  nine  mode  theory  for  laminated  composite  plates.  They  noted  41 
independent  nonzero  stiffness  quantities.  Invariance  of  strain  energy  under  translation  of 
the  reference  surface  limited  the  number  of  stiffness  correction  factors  to  a  maximum 
of  nine  including  two  ”in-plane”  stiffness  factors.  An  alternative  approach  has  been  to 
enforce  equality  of  the  shear  energy  expressions  from  an  approximate  laminate  theory 
and  from  a  stress  field  satisfying  equilibrium.  Chow  [1971]  used  this  method  to  obtain 
the  shear  correction  factor  for  an  orthotropic  laminate  of  symmetric  construction. 
Whitney  [1973b]  extended  it  to  the  case  of  an  asymmetric  orthotropic  laminate.  The 
approximate  stress  field  was  obtained  by  assuming  linear  variation  of  the  in  plane 
stresses  over  the  thickness  of  the  plate  and  using  the  equilibrium  equations  to  solve 
for  the  remaining  components.  Bert  [1984]  used  this  procedure  to  obtain  the  shear 
correction  factor  for  a  laminate  version  of  Levinson's  higher  order  theory.  Reissner 
[1972,  1979]  presented  a  procedure  similar  to  the  one  used  for  homogeneous  plates  to 
obtain  the  transverse  shear  stiffness.  Yang  et  al.  [1966]  also  assumed  the  transverse 
direct  stress  component  to  be  zero.  Whitney  [1970]  assumed  the  integral  of  this  stress 
to  be  zero  over  each  layer  and  applied  his  theory  to  problems  of  cylindrical  bending. 

2.2.2.3  Higher  Order  Theories  [HOT]. 

Whitney  and  Sun  [  1973a]  proposed  a  higher  order  theory  with  quadratic 
polynomial  functions  for  the  in-plane  displacement  components  and  linear  for  the 
transverse  in  order  to  include  the  first  antisymmetric  shear  mode  and  allow  for 
nonzero  normal  strain  in  the  transverse  direction.  A  shear  correction  factoT  was  used. 
In  another  application,  to  laminated  shells,  Whitney  and  Sun  [1974]  assumed  the 
in-plane  displacements  to  be  linear  in  the  transverse  coordinate  and  the  transverse 
displacement  to  be  cubic.  In  both  theories,  a  shear  correction  factor  was  used  to 
enhance  accuracy.  Nelson  and  Lorch  [1974]  employed  quadratic  expressions  for  all  the 
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components  of  displacement.  Lo  et  al.  [  1 977]  used  cubic  functions  for  the  in-plane 
displacement  components  and  quadratic  for  the  transverse  to  have  the  transverse  shear 
contributions  from  the  two  systems  to  be  of  the  same  order  in  the  transverse 
coordinate.  Levinson  [1980]  proposed  using  a  cubic  polynomial  for  the  in-plane 
displacement  components  and  constant  transverse  deflection  for  an  homogeneous,  isotropic 
plate.  To  satisfy  the  stress-free  conditions  on  the  faces  of  the  plate,  the  quadratic  term 
had  to  be  dropped  and  the  coefficient  of  the  cubic  term  was  related  to  that  of  the 
first  order  term  and  the  transverse  component  of  displacement.  Bert  [1984]  extended 
this  idea  to  laminates.  Noting  that  the  equilibrium  equations  in  Levinson's  formulation 
were  inconsistent  in  the  sense  of  a  variational  principle,  Reddy  [1984]  derived  the  field 
equations  from  a  virtual  work  equality.  This  simple  higher  order  theory  involves  the 
same  number  of  variables  as  the  first  order  shear  deformation  theory  but  gives 
parabolic  transverse  shear  strain  variation  through  the  thickness  of  the  plate  and  also 
satisfies  the  stress-free  conditions  at  the  top  and  the  bottom  faces.  This  theory  was 
used  [Bert  1984;  Reddy  1984]  to  analyze  angle-ply  and  cross-ply  laminated  plates  and 
shown  to  be  more  accurate  than  the  first  order  theory.  Table  1  summarizes  some  of 
the  assumed  displacement  fields  used  in  higher  order  theories.  Use  of  higher  order 
terms  in  the  representation  of  displacement  components  increases  the  accuracy  but  the 
additional  terms  introduced  into  the  problem  formulation  increase  the  cost  and 
complexity  of  the  solution  process.  Furthermore,  the  polynomial  representation  over  the 
entire  thickness  of  the  plate  cannot  properly  model  the  discontinuities  in  displacement 
gradients  that  generally  exist  at  the  interfaces  between  layers  with  different  orientation 
of  the  material  axes.  Therefore,  shear  deformable  theories,  whether  FSDT  or  higher 
order,  cannot  take  into  account  the  effects  of  local  deformation  accurately.  This 
difficulty  is  more  pronounced  when  the  shear  rigidities  of  the  constituent  layers  are 
quite  different  [Whitney  and  Sun,  1973;  Lo  et  al.  1977].  Also,  these  cannot  satisfy 
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Table  1:  THEORIES  ALLOWING  FOR  SHEAR  DEFORMATION 


Reference 

In-plane 

Out-of -plane 

Mormal 

Correction 

Force 

Displacement  Displacement 

Strain 

Factors 

Resultants 

Reissner 

Linear 

Constant 

Zero 

1 

8 

[1945] 

4 

1 

Mindlin 

Linear 

Constant 

Zero 

1 

8 

[1951] 

4 

1 

Naghdi 

Linear 

Quadratic 

[1957] 

4 

3 

Whitney 

Quadratic 

Linear 

Constant 

5 

14 

[1973] 

6 

2 

Whitney 

Linear 

Quadratic 

Linear 

[1974] 

4 

3 

Mel son 

Quadratic 

Quadratic 

Linear 

<  9 

17 

[1974] 

6 

3 

Lo 

Cubic 

Quadratic 

Linear 

20 

[1977] 

8 

3 

Levinson 

Cubic 

Constant 

Zero 

1 

8 

[1980] 

4 

1 

Reddy 

Cubic 

Constant 

Zero 

13 

[1984] 

6 

1 
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interlayer  traction  continuity.  None  of  the  existing  higher  order  theories  based  on 
assumed  displacements  can  properly  model  interfacial  stresses,  especially  near  the 
traction-free  edges. 

22.2.4  Discrete  Laminate  Theories  [DLT]. 

In  a  laminated  composite  plate,  the  deformation  depends  upon  the  stacking  sequence 
of  layers  and  upon  the  material  properties  of  each  [Sun  and  Whitney,  1973].  In  order 
to  allow  for  this  properly,  discrete  laminate  theories  in  which  a  layerwise  linear 
variation  of  the  in-plane  displacement  components  is  employed  have  been  proposed  [e.g.. 
Sun  and  Whitney,  1973;  Srinivas,  1973;  Seide,  1980],  Typically,  each  layer  is  treated  as 
an  homogeneous,  anisotropic  plate  and  the  governing  equations  of  all  these  plates  are 
coupled  through  interlaminar  continuity  equations.  Evidently,  to  achieve  greater  accuracy 
of  representation,  it  is  possible  to  divide  each  lamina  into  a  number  of  layers.  Srinivas 
[1973]  assumed  a  piecewise  linear  variation  in  the  in-plane  displacements  with  respect 
to  the  transverse  coordinate  and  the  transverse  displacement  was  assumed  to  be  constant 
over  the  thickness.  This  gave  2N+3  field  equations  for  the  N  layers.  The  theory  did 
not  enforce  continuity  of  traction  across  interfaces.  Sun  and  Whitney  [1973]  used 
kinematic  assumptions  similar  to  Srinivas'.  However,  enforcement  of  continuity  of  shear 
stresses  at  interfaces  reduced  the  number  of  field  variables  to  the  same  as  in  the  first 
order  shear  deformation  theory  regardless  of  the  number  of  layers.  The  degree  of  the 
differential  equation  was  seen  to  increase  with  N,  the  number  of  layers.  Seide  [1980], 
unlike  Srinivas  [1973]  and  Sun  and  Whitney  [1973],  did  not  use  Hamilton's  principle 
but  directly  combined  the  equilibrium  equations  of  the  layers  using  continuities  of 
traction  and  displacement  at  the  interfaces. 

Murakami  [1986]  proposed  a  theory  in  which  layerwise  zigzag  functions  were 
added  to  the  in-plane  displacements  adopted  in  the  ESDT  to  account  for  local  shear 
deformation  effects.  This  theory  was  shown  to  yield  more  accurate  in-plane  stress 
predictions  for  symmetrically  laminated  thick  orthotropic  plates.  Keddy  [1987,  1989] 
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proposed  a  layerwise  theory  as  a  generalization  of  two-dimensional  theories  by  writing 
displacement  components  as  piecewise  smooth  functions  of  the  transverse  coordinate  in 
line  with  finite  element  concepts. 

2.2.3  Stress-based  Theories. 

Pagano  [l977a,b;  1978a, b]  proposed  a  theory  which  assumed  linear  variation,  with 
respect  to  the  transverse  coordinate,  of  the  in-plane  stresses.  The  remaining  stress 

components  were  derived  from  the  three-dimensional  equilibrium  equations.  A 
variational  principle  was  used  to  develop  constitutive  equations  relating  force  resultants 
to  ’’generalized  displacements.”  Continuity  of  tractions  as  well  as  displacements  at 

interfaces  was  explicitly  satisfied.  The  boundary  value  problem  was  stated  in  terms  of 
the  generalized  displacements  and  interfacial  tractions.  The  only  approximation  made  in 
this  theory  is  the  assumption  of  linear  variation  of  the  in-plane  stress  components  over 
each  layer  thickness.  As  each  layer  can  be  subdivided  into  as  many  sublayers  as  one 
wishes,  the  theory  ought  to  lead  to  a  sequence  of  solutions  of  the  problem  which 
would  converge  to  the  exact  solution  with  reduction  in  the  sublayer  thickness. 

In  Pagano's  theory,  the  number  of  simultaneous  partial  differential  equations  to  be 
solved  is  13N.  This  has  been  the  principle  difficulty  in  applying  this  theory  to 

practical  problems.  To  overcome  this,  Pagano  and  Soni  [1983]  proposed  a  ’’global-local" 
model  in  which  a  laminate  would  be  divided  into  a  local  zone  and  a  global  portion.  A 
higher  order  displacement-based  theory  would  be  used  to  solve  the  global  problem  and 
then  Pagano's  theory  would  be  applied  to  the  local  zone  for  better  accuracy  in  that 
region. 

Green  and  Naghdi  [1981,  1982]  proposed  a  generalized  discrete  laminate  theory 

incorporating  thermal  and  dynamical  effects.  Table  2  gives  a  comparison  of  the 
discrete  laminate  theories. 
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Table  2:  DISCRETE  LAMINATE  THEORIES 


Reference 

Type 

In-plane  Transverse  Field 

Displacements  Displacement  Variables 

Seide 

[1980] 

Displacement 

Linear 

Constant 

2N  ♦  3 

Sun  and  Whitney  Displacement 
[1973] 

Linear 

Constant 

2H  ♦  3 

Srinivas 

[1973] 

Displacement 

Linear 

Constant 

211  ♦  3 

Pagano 

[1977,1978] 

Stress 

Linear 

Quadratic 

Shear 

Cubic  Hormal 
Stress 

13N 

13 


2.2.4  Free-Edge  Delamination  Specimens. 

Free-edge  delamination  specimens  are  tested  under  uniform  axial  strain.  This  implies 
that  the  axial  component  of  displacement  is  linear  in  the  longitudinal  coordinate.  Away 
from  the  grips,  it  is  reasonable  to  assume  that  the  other  two  displacement  components 
are  independent  of  the  longitudinal  coordinate.  This  reduces  the  three-dimensional 
problem  to  one  involving  three  variables  (components  of  the  displacement  vector)  but 
these  are  functions  of  only  two  independent  spatial  variables,  viz.,  the  transverse 
coordinate  and  the  ’’other”  of  the  in-plane  coordinates. 

Specialization  of  the  plate  theories,  described  earlier,  to  the  case  of  a  free-edge 
delamination  specimen  reduces  the  dimension  of  the  problem  from  two  to  one.  Thus, 
the  dependence  of  various  field  variables  is  on  the  transverse  in-plane  coordinate  only. 

2.3  FINITE  ELEMENT  MODELLING. 

Even  highly  simplified  models  of  laminated  composites  lead  to  very  complex  sets 
of  equations.  Exact  solutions  to  these  equations  have  been  passible  only  for  cases  of 
simplest  geometry  and  load  configurations.  Use  of  numerical  methods  to  obtain 
approximate  solutions  to  the  field  equations  is  thus  imperative.  Both  finite  difference 
and  finite  element  methods  have  been  employed. 

Hulbert  and  Rybicki  [1971]  used  a  boundary  point  least  squares  technique.  Pipes 
and  Paganao  [l 970]  used  the  finite  difference  method.  Most  of  the  finite  element 
studies  have  been  for  nlane  stress  or  three-dimensional  models  of  composite  laminates. 
Rybicki  [1971]  used  Maxwell's  stress  functions  approximated  as  fifth  degree  Hermite 
polynomials  to  get  a  stress  field  satisfying  equilibrium.  This  gave  648  unknowns  per 
element.  Minimization  of  the  complementary  energy  was  used  to  evaluate  the 
'■oefficients  for  the  stress  functions.  Isakson  and  Levy  [1971]  used  a  shear  interlayer 
element  and  a  displacement-based  formulation.  Constant  strain  triangles  were  used  to 
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model  the  laminate  layers.  Wang  and  Crossman  [1977a,  1977b]  used  constant  strain 
elements  and  studied  a  plane  strain  slice  of  an  angle-ply  laminate  specimen.  Rybicki 
and  Schmueser  [l 978]  used  Wilson's  SAP  IV  program  to  carry  out  a  three-dimensional 
analysis  based  on  eight-point  ’’brick”  isoparametric  elements.  Most  investigators 
confirmed  Puppo  and  Evensen's  [1970]  and  Pagano  and  Pipes  [1973]  findings.  However, 
Bar-Yoseph  [1983]  noted  that  Whitcomb  et  al.'s  [1982]  findings  are  in  contradiction  of 
Pipes  and  Pagano's  [1970]  finite  difference  solution  and  the  displacement-formulated 
perturbation  solution  by  Hsu  and  Herakovich  [1977].  All  of  the  analyses  were  based  on 
linear  elastic  material  behavior  and  constant  axial  strain.  None  could  allow  for 
cumulative  damage.  Nishioka  and  Atluri  [1980,  1982]  used  an  assumed  stress  approach. 
Spilker  and  Chou  [1980a, b]  presented  a  special  purpose  hybrid-stress  multi-layer  finite 
element  satisfying  traction-free  edge  conditions  exactly  and  applied  it  to  a  symmetric 
cross-ply  laminate.  Mixed  (displacement-stress)  formulations  have  been  attempted  [e.g., 
Labbe  et  al,  1982], 

To  admit  crack-tip  singularity  into  the  analysis,  special  singularity  elements  have 
been  used  for  some  time.  Wang  et  al.  [1975a,  1975b],  and  Wang  and  Mandell  [1977] 
used  a  two-dimensional  hybrid  stress  element  including  a  crack-tip  singularity  element 
embedded  in  a  matrix  interlayer  between  plies  of  the  laminate.  Muskhelishvilli's 
solution  was  the  basis  for  the  special  element.  Accuracy  of  one  percent  for  the  stress 
intensity  factor  calculation  wc  claimed.  L’eng  et  al.  [  1 977]  applied  this  procedure  to 
study  failure  of  a  notched,  unidirectional,  boron-epoxy  composite  under  simple  tension 
along  the  fibers  and  perpendicular  to  the  notch.  Raju  et  al.  [1980]  and  Whitcomb  et 
al.  [1982]  tested  the  capability  of  the  finite  element  method  to  solve  for  interlaminar 
stress  in  composite  laminates.  It  was  felt  that  with  sufficient  refinement,  the  method 
could  give  a  highly  accurate  estimate  of  the  interlayer  stress  field  except  in  the  two 
elements  closest  to  a  singularity  or  discontinuity.  These  conclusions  were  based  on  use 
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of  8-point  ’’brick”  (trilinear  lagrangian)  three-dimensional  elements  and  8-point 
isoparametric  quadrilateral  plane  strain  elements.  Stanton  and  Crain  [1980]  sought  to 
economize  on  the  computer  effort.  Linear  constraints  were  applied  to  tricubic 

Lagrangian  three-dimensional  elements  to  get  a  lower  order  displacement  along  selected 

axes  while  retaining  the  higher  order  variation  in  material  properties  and  higher  order 

geometry.  Transition  elements  were  used  to  pass  from  discrete  ply  modelling  to 

composite  laminate  modelling.  Raju  and  Crews  [1981]  used  a  polar  mesh  near  the 
intersection  of  the  interface  and  the  free  edge  along  with  a  log-linear  radial  distance 
from  the  singular  point.  Bar-Yoseph  [1982]  presented  a  method  based  on  composite 
expansion  and  assumed  stress  approach.  Wang  and  Slomiana  [1982]  introduced  finite 
elements  allowing  for  fracture  based  on  Griffith's  criterion  to  simulate  crack  growth. 
This  is  an  extension  of  earlier  work  by  Wang  and  Crossman  [1977a].  All  three  modes 
of  crackin'1  were  considered.  Rybicki  et  al.  [1977]  used  an  energy  release  rate  criterion 
in  finite  element  analysis  of  crack  growth.  Table  3  lists  a  comparison  of  various 
methods  used  for  numerical  solution  of  the  free-edge  delamination  specimen  problem. 

In  order  to  set  up  finite  element  solution  procedures,  variational  formulations  are 
generally  used.  For  analysis  of  laminates,  minimization  of  potential  energy, 
minimization  of  complementary  energy,  Reissner's  mixed  variational  principle,  hvbrid 
variational  formulations,  global-local  formulations  etc.,  have  been  used. 

All  the  investigators  report  success  with  whatever  method  they  adopted.  Good 
agreement  with  test  data  is  claimed.  However,  it  is  noteworthy  that,  in  general,  no 
mention  is  made  of  computational  efficiency.  It  is  well  known  that  triangular  elements 
do  not  give  good  distribution  of  stress  except  in  the  simplest  cases.  Again,  a  small 
number  of  trilinear  elements  will  not  give  good  stress  distribution  which  could  be  used 
as  the  basis  for  simulation  of  initiation  and  progresive  growth  of  cracks  and/or 
delamination.  It  has  also  been  noted  [e.g.,  (’hang  1981;  Jenq  1982]  that  mixed 
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Table  3:  STRESS-DISTRIBUTION  IN  A  FPEE-EDGE  DELAMINATION  SPECIMEN. 


Reference 

Method  of  Analysis 

Stress  Components 

Puppo  k  Evensen 
[1970] 

Approximate  Elastic  Solution 

cr  cr  c r  ,<7  . 

XX,  yy,  xz’  xy 

Pipes  k  Pagano 
[1974] 

Approximate  Elastic  Solution 

cr  ,cr  ,  cr  ,<r 

xx’  xy’  xz’  yz 

Hsu  k  Herakovich 

Perturbation  Technique 

cr  ,cr  ,<r 

zz  xz’  yz 

[1977] 

Pagano 

[1978] 

Stress-Based 

Theory  of  Plates 

cr  y<j  ,cr  ,<r  ,cr  ,cr 

xx’  yy  Z2  yz’  zx’  xy 

Pagano  &  Soni 
[1983] 

Global-Local  Model 

a  ,c r  yCr  ,cr  ,cr  ,cr 

xx  yy  zz  yz  zx  xy 

Pipes  k  Pagano 
[1970] 

Finite  Difference 

cr  ,a  ,c r  ,cr  ,cr  ,c r 

xx  yy  zz’  yz  zx’  xy 

Wang  k  Crossman 

Finite  Element 

cr  ,cr  ,cr  ,c r  ,cr  ,cr 

xx'  yy’  zz’  yz’  zx’  xy 

[1977a] 

Constant  Strain  Triangle 

Finite  Element 

Whitcomb 

8-Noded  Isoparametric 

cr  .cr  ,  cr  , cr  ,cr  .  cr 

xx’  yy’  zz*  yz'  zx'  xy 

[1982] 

Element 

Rybicki 

[1971] 

Finite  Element  Method 
Equilibrium  Stress  Method 

cr  ,  cr  ,  cr  ,  cr  ,  cr  ,  cr 

xx’  yy’  zz’  yz’  zx  xy 

Spilker 

[1980b] 

Finite  Element  Method 

Hybrid  Stress  Model 

cr  ,  cr  ,  cr  ,  cr  ,  cr  ,  cr 

xx  yy’  zz  yz’  zx  xy 

Wang  k  Yuan 
[1983] 

Finite  Element 

Singular  Hybrid  Element 

cr  y  cr  y  cr  t  cr,  ,  o*  ,  cr 

xx  yy  zz  yz  zx  xy 
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(displacement-stress)  formulations  result  in  poorly  conditioned  matrices  and  the  solutions 
have  spatially  oscillatory  errors.  Most  finite  element  formulations  based  on  use  of 
special  singularity  elements  can  adequately  predict  onset  of  fracture  in  the  vicinity  of 

the  crack-tip  but  are  not  convenient  for  prediction  of  crack  propagation  and  arrest. 

Also,  they  cannot  predict  the  initiation  of  fracture  in  an  initially  uncracked  material. 

Some  work  [e.g.,  Pryor  and  Barker  1970;  .Vlawenya  and  Davies  1973;  Hinton  1976; 
Reddy  1982]  has  been  done  on  flexure  of  composite  laminates.  In  general,  the  finite 
element  interpolation  along  the  transverse  axis  is  chosen  to  follow  one  of  the  plate 
theories.  In-plane  variation  follows  the  usual  finite  element  procedures.  Lagrangian  as 
well  as  isoparametric  interpolations  have  been  used.  Kngblom  and  Ochoa  [1985]  used 
through-the-thickness  elements  for  restricted  applications.  Mawenya  and  Davies  [1973] 
used  a  layerwise  model.  Recently,  Reddy  [1982,  1987]  and  Kuppuswamy  and  Reddy 

[1984]  have  used  the  layerwise  theory  for  laminated  plates. 

2.4  THE  PRESENT  RESEARCH. 

Review  of  the  existing  literature  on  the  subject  showed  that  a  proper  theory  for 
modelling  bending  and  stretching  behavior  of  composite  laminate  plates,  with  the 
capability  of  providing  reliable  solutions  for  interlaminar  stresses,  was  lacking.  The 
available  theories  were  based  or:  a  number  of  ad-hoc  assumptions  and  were  applicable 
to  only  certain  special  situations.  Pipes  and  Pagano  [1970]  pointed  out  that  the 
interlaminar  normal  stress  is  primarily  responsible  for  delamination  failures.  It  was 
extremely  important  to  develop  adequate  theoretical  models  which  could  be  used  to 
predict  stresses  at  the  interfaces  with  the  required  degree  of  accuracy.  Finite  element 
models  ensuring  continuity  of  tractions  across  material  interfaces  are  rare.  Fot  plate 
theories  there  was  no  concensus  regarding  consideration  of  shear  deformation  and  most 
theories  did  not  satisfy  interlayer  continuity  of  tractions  and  were  inadequate  for  study 
of  local  effects  causing  damage.  Pagano’s  theory,  based  on  assumption  of  linear 
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variation  of  the  in-plane  stress  components  over  each  layer,  has  given  excellent  results. 
However,  its  application  has  been  restricted  to  a  few  layers.  It  appeared  necessary  to 
develop  procedures  w'hich  would  permit  application  of  Pagano's  theory  to  multi-layer 
laminate  configurations.  This  led  to  the  following  lines  of  investigation  for 
development  of  the  theoretical  models: 

1.  Systematic  development  of  models  of  combined  bending  and  stretching  of 
plates  properly  allowing  for  interlayer  continuity  of  tractions,  including  revi¬ 
siting  existing  theories; 

2.  Finite  element  implementation  of  the  theoretical  models  developed  under  (l) 
above; 

3.  Development  of  finite  element  models  for  analysis  of  free-edge  delamination 
specimens. 

Under  item  1  above,  an  hierarchical  approach  to  the  construction  of  assumed 
displacement  type  theories  of  laminated  plates  was  identified.  The  constitutive 
equations  in  each  case  were  derived  from  a  general  variational  principle.  These  theories 
do  not  require  the  use  of  ad-hoc  correction  factors  but  directly  and  explicitly  contain 
the  influence  of  the  layup  and  the  layer  properties  upon  the  constitutive  relations. 
Pagano's  theory  was  restated  in  a  form  convenient  for  application  of  finite  element 
procedures.  It  was  seen  that  the  constitutive  equations  for  this  theory  need  not  be 
based  upon  any  variational  principle  but  could  be  derived  directly  by  integrating  the 
material  constitutive  relations.  The  first  order  coupled  shear  deformation  theory  and  the 
restated  Pagano's  theory  were  implemented  in  finite  element  computer  programs. 
Additionally,  approximation  of  a  free-edge  delamination  specimen  as  a  segment  of  an 
annulus  with  a  large  radius  was  implemented  in  a  finite  element  computer  program.  A 
significant  result  was  the  development  of  continuous  traction  finite  elements  using  a 
cubic  displacement  representation  for  analysis  of  free-edge  delamination  specimens.  In 
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each  case  the  influence  of  element  size  (mesh  refinement)  on  accuracy  was  carefully 
examined.  Wherever  possible,  microcomputer  versions  of  the  codes  were  prepared.  All 
the  computer  codes  were  carefully  documented.  The  codes  have  a  modular  structure  to 
permit  easy  update  and  enhancement  in  the  future.  The  validity  of  the  proposed 
models  was  checked  against  available  exact  solutions  and  test  data.  The  computer  codes 
were  used  to  predict  stresses  and  deformations  in  multi-layer  free-edge  delamination 
specimens. 
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SECTION  III 


EQUATIONS  OF  THREE-DIMENSIONAL  LINEAR 
ELASTICITY  AND  THEIR  SPECIALIZATION  TO 
FREE-EDGE  DELAMINATION  SPECIMENS 

3.1  INTRODUCTION. 

We  start  by  summarizing  the  well-known  equations  of  three-dimensional  linear 
theory  of  elasticity.  We  shall  then  specialize  these  for  the  case  of  free-edge 
delamination  specimens  allowing  for  the  special  geometry  and  loading  configuration. 
These  consist  of  the  fact  that  the  specimens  are  uniformly  stretched  in  the  longtudinal 
direction  and  that  the  laminate  possesses  certain  symmetries. 


3.2  THREE-DIMENSIONAL  ELASTICITY. 
3.2.1  Kinematics. 


If,  using  a  rectangular  cartesian  frame  of  reference,  u(  are  the  components  of  the 
vector  describing  the  displacement  of  an  arbitrary  particle  in  an  elastic  solid,  the 
cartesian  components  eu  of  the  infinitesimal  strain  tensor  are  given  by 


e  =  -L  (u  +  u  ) 
>j  2  *-J 


(1) 


Here,  and  in  the  sequel,  we  use  standard  index  notation.  Roman  indices  take  on  the 
range  of  values  from  1  to  3.  Summation  on  repeated  indices  is  implied  unless 
otherwise  indicated.  A  subscripted  comma  denotes  differentiation  with  respect  to  the 
spatial  coordinate  indicated  by  the  subscript  following  the  comma.  Parentheses  around  a 
pair  of  indices  denote  the  symmetric  prrt  of  the  quantity  with  respect  to  the  two 
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indices.  Greek  indices  take  on  values  1  and  2.  Parentheses  around  a  single  index 
denote  "no  sum”  on  that  index. 


The  strain  tensor,  defined  by  (1)  is  symmetric,  i.e. 

e  =  e  =u,  , 

>j  p  M 


(2) 


3.2.2  Equilibrium. 

If  f,  are  the  components  of  the  body  force  vector  acting  upon  an  arbitrary  unit 
volume  of  the  solid,  the  static  equilibrium  equations  for  the  solid  are: 

cr  +  f  =  0  over  R  (3) 

ij.'  j 

Here  c ri}  are  the  components  of  the  symmetric  Cauchy  stress  tensor,  R  is  the  deformed 
configuration  of  the  solid  on  which  equilibrium  is  defined.  However,  for  linear  theory, 
no  distinction  is  made  between  the  undeformed  and  the  deformed  state  for  the  purpose 
of  defining  the  stress  components  and  the  equations  of  equilibrium.  Often,  the  quantity 
on  the  left-hand  side  of  (3)  is  referred  to  as  the  "disequilibrium  force  vector.”  Using  d, 
to  represent  components  of  the  disequilibrium  force  vector,  their  vanishing  everywhere 
implies  pointwise  equilibrium.  However,  if  we  consider  the  inner  product 


P  = 


(4) 


the  product  vanishes  for  every  bounded  function  gt  if  and  only  if  the  quantity  d,  is 
zero  almost  everywhere.  Further  noting  that  any  bounded  function  can  be 
approximated  as  closely  as  desired  by  a  polynomial  with  rational  coefficients,  vanishing 
of  d;  can  be  replaced  by  P  =  0  for  g,,  a  polynomial  with  rational  coefficients  and  of 
arbitrary  degree  in  the  spatial  coordinates.  This  weak  form  is  the  basis  of  theories  of 
mechanical  behavior  of  surface  structures  e.g.,  plates  and  shells.  We  shall  refer  to  this 
concept  again  in  Section  IV. 
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3.2.3  Constitutive  Relations. 


For  a  linear  elastic  material,  the  constitutive  relations  have  the  form 


ffij  ^Mijekl 


(5) 


Here  E1Jki  are  components  of  the  isothermal  elasticity  tensor.  Because  of  symmetry  of 
the  strain  tensor, 

E.  ,  =  E ....  (6) 

gkl  jikl 

If  there  are  no  body  force  couples,  the  stress  tensor  is  symmetric.  This  leads  to 

E  v,  =  E  .  (7) 

ijkl  ijlk 


Further,  if  it  is  assumed  that  the  stress  is  independent  of  the  strain  path  and  an 
energy  density  function  of  strain  exists  such  that  each  component  of  the  stress  tensor 
is  the  gradient  of  this  function  along  the  ’’corresponding”  component  of  the  strain 
tensor. 


(8) 


Thus,  the  number  of  independent  constants  defining  the  elastic  material  reduces  to  21. 
For  a  material  having  symmetry  about  a  plane  the  number  reduces  to  13  and,  for  an 
orthotropic  material,  the  number  is  9.  The  relationship  can  then  be  written  in  the 
form 


cr  =  E  ,  „e  ,  +  E,,  „e,, 

a/3  yoa/3  yo  33a0  33 

(9) 

a  .  "2E  ,  ,e  , 

a3  y3a3  >3 

(10) 

cr  =  F  e  +  E  e 

33  y633  y6  ~  ^3333  C33 

(11) 

inverse  relationships  are 

eo0  =  Cyf.o/3  ^yb  +  ^'op 33  ^33 

(12) 

e  =  2C  ,  , or 

o3  >3or3  >3 

(13) 

e,,=C  , ,,  cr  ,  +  C  cr 

33  y>03  yb  3333  33 

(14) 
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3.3  SPECIALIZATION  TO  FREE-EDGE  DELAMINATION 


SPECIMENS. 


3.3.1  The  Free-Edge  Delamination  Specimen. 

Figure  2  shows  a  symmetric  laminate  specimen  under  a  uniform  axial  strain  e0  . 
For  this  case,  away  from  the  ends,  the  displacement  field  for  any  transverse  section 
(x,  =  constant)  has  the  following  form  [Pipes  and  Pagano  1970] 
u|(xi)  =  e0x,  +U1(x,,x3) 

u,(x()  =  U,(x2,x3)  (16) 


and 


u3(xt)  =  U3(x2,  x3) 


(17) 


3.3.2  Effect  of  Specimen  Symmetries. 

Symmetry  of  loading  as  well  as  geometry  about  the  midplane  (x3  =  0)  and 
rotational  symmetry  about  the  axis  x3  through  the  middle  point  result  in  the  following 
relationships: 


Ua(x2,  x3)  Uo(x2,x3) 

U3(x2 ,  -  x3)  =  -  U3(x2 ,  x3) 

L /  ~  X0  ’  X3^  —  ~  ^  ’  X3^ 

U3(  -x2,x3)  =  U3(x2,x3) 


(18) 

(19) 

(20) 
(21) 


Using  the  chain  rule  of  differentiation,  (18),  (19)  and  (21)  yield 
-  Uq3(x2  ,  -  x3)  =  Ua3(x2 ,  x3) 

^3,2^X2  ’  _  X3^  ^3,2^X2  ’  X3^ 

—  U'3 2(  —  x2, x3)  -  UJ2(x2,x3) 


(22) 

(23) 

(24) 


Setting  x,  =  0  in  (19),  (22),  and  (23), 
U3(x2,0)  =  0 


(25) 
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(26) 


U 3  2(x , ,  0)  =  Uq  3(x  2,0)  =  0 
Equation  (26)  also  gives 

(U3.2-Ua.3)(x2.0)  =  0  f0r  a11  *2 

Equation  (20)  gives 

b’2(0,x3)  =  0  for  all  x3 


(27) 

(28) 


and,  as  a  consequence. 


U23(0,x3)  =  0  for  all  x3 

(29) 

Equation  (24)  gives 

L32(0,x1)  =  0  for  all  x3 

(30) 

Equation  (29)  along  with  (30)  leads  to 

-y23  =  0  for  all  x3  at  x2  =  0 

(31) 

and 

<u«-uuW’° 

(32) 

Summarizing,  for  the  free-edge  delamination  specimen, 

the  following  constraints 

apply: 

U2(0,x3)  =  U3(x2,0)  =  0 

(33) 

^  3,2  ~  ^2.3^(0.x3)  “  ^  3.2  ~  ^or^XyO)  ~  ® 

(34) 

y23U2,0)  =  y2i(0,\3)  =  0 

(35) 
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SECTION  IV 


EQUATIONS  GOVERNING  BENDING  AND 
STRETCHING  OF  PLATES 

4.1  INTRODUCTION. 

To  reduce  the  three-dimensional  problem  of  a  plate  to  one  in  two  dimensions,  the 
dependence  of  the  field  variables  on  the  transverse  coordinate  is  made  explicit.  This 
involves  assumptions  regarding  dependence  of  stress  or  displacement  on  the  coordinate 

X3‘ 

In  the  present  research  effort,  using  the  discrete  laminate  concept,  two  approaches 
were  considered.  In  one,  assumptions  were  made  regarding  displacement  variation  with 
respect  to  the  coordinate  x3.  An  innovative  component  of  the  present  effort  was  that 
the  constitutive  relations  were  derived  directly  from  a  variational  theorem  using  a 
generalization  of  Reissner's  work  [1984].  This  revealed  the  existence  of  a  coupling 
between  the  force  resultants  and  deformations  in  various  layers. 

In  the  second  approach,  originally  proposed  by  Pagano  [1978],  assumptions  were 
made  regarding  stresses  which  were  required  to  be  in  pointwise  equilibrium.  The 
stress-strain  relations  were  satisfied  pointwise  as  well.  The  strains  were  integrated  using 
appropriate  weighting  functions  to  yield  constitutive  relations  for  the  force  resultants 
in  terms  of  generalized  displacements  and  their  gradients.  Continuity  of  displacements 
as  well  as  of  tractions  was  enforced  across  interlayer  surfaces.  No  assumptions  were 
necessary  regarding  dependence  of  the  displacement  components  upon  x3. 

In  this  Section  we  summarize  the  first  approach.  Details  are  given  in  reports  listed 
as  items  B.1.1;  B.1.2  and  dissertations  listed  as  items  B.4.2  and  B.4.5  in  Appendix  B. 
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The  second  approach  is  discussed  in  Section  V. 


4.2  GOVERNING  EQUATIONS. 


4.2.1  Kinematics. 

In  general,  using  a  polynomial  representation  of  the  dependence  of  the  components 
of  displacement  upon  the  coordinate  x3  through  the  thickness  of  a  plate,  one  could 
write 

uXXj)  =  wfxp)  +  x3  <f>.(xp)  +  *,(xp)  +  ....  (36) 

Most  theories  of  bending  and  stretching  of  plates  employ  a  polynomial  up  to  a  certain 
degree.  For  example,  using  terms  up  to  the  first  order  for  the  in-plane  displacement 
components  and  only  the  constant  term  for  u3  gives  the  first  order  theory  described  by 
the  assumptions 

(37) 

(38) 

Corresponding  to  these  displacements,  the  strains  aie: 

(39) 

y  .  =  v3  +<(>  (40) 

•  a  3  3, a 


Uc(Xi)  =  VW  +  X3W 
U3(X,)  =  V3(V 


e^  =  V(^)(x>)  +  X3WX>) 


e33  =  0 


(41) 


For  this  case,  the  quantity  <f>laJ})  is  referred  to  as  the  curvature  of  the  plate.  Using 
the  general  polynomial  representation,  the  components  of  strain  are 

ea/5  =  V(aJXy)  +  X3  +  X3  + -  (42) 

ya3  =  v3  o  +  x303a  +  x3 1 l/3a  +  X3  XXo  + - +  +  2  x3  \l/a  +  3  X3  xa  +  .  (43) 

e33  =  V3,3  +  2  X3  ^3  +  3  X3  *3  +  —  (44) 

For  the  discrete  laminate  theory,  the  above  equations  apply  to  each  layer. 
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4.2.2  Equilibrium. 

In  line  with  the  concepts  outlined  in  Section  4.1,  the  ’’disequilibrium  force  vector” 
in  Section  3.2.2  is  orthogonal  to  an  arbitrary  function  over  the  thickness  of  the  layer. 
Noting  that  any  bounded  function  over  a  finite  domain  can  be  approximated  by  a 
polynomial  of  sufficiently  high  degree  and  with  rational  coefficients,  the  condition  for 
equilibrium,  using  a  polynomial  in  x3  is 


+  fpx"dx3  =  0,  for  n  =  0,  1,  2 . 


(45) 


I7or  a  cylindrical  solid  with  area  A  in  coordinates  xa  and  thickness  t  in  coordinate  x3, 
(45)  is 


t 


2 

/  / (tr,>,  +  fJ)X3dX3dA=0 

A 

2 

Here  t  could  be  a  function  of  xa  If  we  require 
(46)  vanish  over  arbitrary  regions  in  the  domain 


(46) 

that  the  quantity  on  the  left  side  of 
defined  by  x„ 


/(cr  +f)x"dx  =0 
v  i>i  /  3  3 


For  n  =  0,  this  leads  to  the  equations 


N'  +  cr  —  cr  +  F  „  =  0 

Ckfi.a  /33  P  3  p 


0  +0-  — cr  +F,  =  0 

33  33  3 

Here  we  have  introduced  force  resultants 


(47) 


(48) 

(49) 


(cr  „,cr  „f  )dx 

a/3  a3  i 


3 


2 


(50) 
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The  superscripted  +  and  -  ,  respectively,  denote  the  top  and  the  bottom  surfaces  of  the 
plate.  Considering  n  =  1  and  only  the  equilibrium  equation  corresponding  to  i  =  1,2 
we  get 

where  we  define  additional  generalized  forces 

t 

2 

3^X3  ^2) 

_  t 
2 

The  quantities  Q»,  NuS,  Mtta,  are,  respectively,  the  shear  force,  the  normal  force  and  the 
bending  moment  resultants  for  the  plate. 

We  note  that  for  a  consistent  theory,  the  number  of  equilibrium  equations  must 
equal  the  number  of  free  kinematic  field  variables.  Thus,  for  the  first  order  theory 
represented  by  the  kinematic  assumptions  of  (37)  and  (38),  five  equations  of 
equilibrium  are  required.  Equations  (48)  ,  (49)  along  with  (51)  provide  these.  If  ua  is 
quadratic  and  u3  linear  in  x3,  there  are  eight  kinematic  variables  requiring  eight 

equations  of  equilibrium.  Following  the  argument  given  above,  we  shall  use  i  =  1,  2, 
3  for  n  =  0,  1  along  with  i  =  1,  2  for  n  =  2  to  furnish  the  eight  equations. 

Extending  the  argument  to  the  case  of  uQ  cubic  and  u3  quadratic  in  x3,  we  have 

eleven  kinematic  variables  and  the  corresponding  equations  of  equilibrium  will  be  set 
up  using  i  =  1,  2,  3  for  n  =  0,  1,  2  and  i  =  1,  2  for  n  =  3.  This  provides  a 

systematic  approach  to  an  hierarchical  development  of  higher  order  theories  of  any 
order  desired. 
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4.3  CONSTITUTIVE  EQUATIONS  FOR  ELASTIC  PLATES. 


4.3.1  Introduction. 

In  setting  up  equations  of  equilibrium,  force  resultants  i.e.,  quantities  of  the  type 

X_ 

2 

f  <V3dX3  (53) 

_-L 

2 

were  introduced.  For  n  =  0,  there  are  five  quantities  viz.,  N0(S,Qtt.  For  n  =  1,  there  are 
six  additional  quantities,  viz.,  Mafi,  corresponding  to  i  =  1,  2  and  Ma},  N33  defined  by 


Mo3  =  /°'o3X3dX3  ^54) 

~  2 


2 

N3J  =  /°-l3dX3 

-JL 

2 

For  n  =  2  and  i  =  1,  2  there  would  be  additional  quantities 

2 

Pa0=/%SX  3dX3  (56) 

__t 

2 

The  number  of  force  resultants  is  invariably  larger  than  the  number  of  equations  of 
equilibrium.  For  this  reason,  the  equations  of  equilibrium  are  not  sufficient  to  define 
the  resultants.  However,  if  the  equations  of  equilibrium  are  expressed  in  terms  of  the 
kinematic  variables,  the  number  of  equations  of  equilibrium  being  equal  to  the  number 
of  kinematic  variables,  a  solution  could  exist.  The  relationships,  expressing  force 
resultants  in  terms  of  kinematic  variables,  are  the  constitutive  equations  for  the 
problem.  For  the  first  order  theory,  for  five  equations  of  equilibrium,  there  are  eight 
force  resultants  needing  as  many  constitutive  relationships.  For  uo  quadratic  and  u3 
linear  in  x3  there  are  eight  equations  of  equilibrium  in  14  force  resultants  needing  14 
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constitutive  relations,  and  so  on. 

For  a  laminate  there  are  kinematic  variables,  equilibrium  equations  and  force 
resultants  associated  with  each  layer.  In  fact,  additional  field  variables  appear  in  the 
form  of  interlayer  tractions  which  are  a  priori  unknown  and  are  of  special  importance 
in  their  influence  on  damage  mechanisms.  Corresponding  to  these  interlayer  tractions, 
equations  of  balance  arise  as  continuity  conditions  for  the  components  of  displacement. 
For  preexisting  and  specified  delamination  opening  or  slip,  these  would  be  the  specified 
measures  of  discontinuity.  Thus,  the  number  of  kinematic  variables  for  the  first  order 
theory  is  5N,  where  N  is  the  number  of  layers,  and  there  are  8N  force  resultants 
along  with  3(N  -  l)  interlayer  tractions  with  as  many  displacement  continuity 
conditions.  For  a  higher  order  theory  based  on  uQ  quadratic  and  u3  linear  in  x3,  there 
are  8N  kinematic  variables  with  as  many  equilibrium  equations,  14N  force  resultants 
along  with  3(N  -  1)  interlayer  tractions  and  as  many  displacement  continuity  equations. 

To  derive  these  constitutive  relationships,  one  procedure  would  be  to  simply 
integrate  the  stress-strain  relations  for  the  material.  Another  is  to  proceed  via 
appropriate  variational  formulations  written  for  three-dimensional  elasticity  and  then 
specialized  to  reflect  the  kinematic  assumptions  and  the  definitions  of  the  force 
resultants. 

Herein,  we  describe  the  traditional  approach.  An  alternative,  suggested  by  Reissner's 
work,  used  for  the  development  of  a  hierarchical  approach  to  constitutive  theory  of 
plates  is  discussed  in  Section  V. 


4.3.2  Direct  Integration  of  Material  Constitutive  Relations. 

The  general  linear  elasticity  relationships  described  by  (5)  upon  specialization  to  the 
case  of  a  lamina  with  material  symmetry  about  the  middle  surface,  reduce  to  (9) 
through  (11)  which,  for  the  kth  lamina,  are: 


(k)  __  t/k)  (k)  ,  j,(k)  (k) 

%  ~  tyl  +  t’33or^  e33 


(57) 
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(58) 


00  «  ,4k)  (k) 

<rc.3  =  2Ey3Q3ey3 


_  ,4k)  (k)  ,  ,4k)  (k) 

^33  ~  Ey633  ey6  +  E3333  e33 


(59) 


In  developing  theories  of  thin  homogeneous  plates,  it  is  customary  to  eliminate  e33 
between  (57)  and  (59)  and  then  to  assume  a  condition  of  plane  stress  to  get 

_(k)  _  ,4k)  _(k)  (60) 


_  p<k) 

^  ybofi^yb 


where 


p(k)  p(k) 

p<k)  _  p(k)  ty?,33r'33o^ 

y&a/3  y  &o/3  (k) 

E3333 


(61) 


Noting  (42)  and  (43),  the  material  constitutive  relations  upon  appropriate  integration  for 
the  first  order  theory  and  an  homogeneous  layer  give  the  following  constitutive 
relationships  for  the  force  resultants: 

(62) 


N(k)  =  t  E<k)  v(k) 

lk  y6a0  V(y.6) 


...ik).  *k  c(k)  ,(k) 

Map  12  y6o,0  ^(y.6) 

Q!k)=\EyL(<Xk)) 


(63) 


(64) 


We  note  that  the  elimination  of  e33  is  not  consistent  with  the  kinematic  assumption  of 
vanishing  e33.  For  the  higher  order  theories  in  which  e33  is  nonvanishing,  the 
assumption  of  generalized  plane  stress  could  be  used.  However,  for  laminates,  the 
transverse  stress  at  the  interface  is  not  negligible  and  indeed  governs  delamination.  Here 
this  assumption  would  not  be  admissible.  The  general  relationships  based  on  (42) 
through  (44)  would  be 


N(k)  =  t  E<k) 

lk  y&ag 


<k) 


lk  .i<k> 


,  ,  K  |  IK/  , 

V(y.6)  +  72  ^(>'6)  +  ' 


+  t  E<k) 

T  lk  C33aP 


(k)  ^  \  Jk) 


M 


v  +  —  v  + 

V 3,3  4  *3 

_,3c(k)  (  1  ,(k),  7  ,  ,3,4k)  (  1  ,(k),  A 
=  lkEy6o  A  TI  . )  lkE33a  A  X  *3  + ) 


(65) 

(66) 
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etc.  The  assumption  of  generalized  plane  stress  was  first  used  by  Mindlin  [l95l]  for 
homogeneous  plate  theory  and  later  by  Whitney  {1970]  for  laminated  plates.  Wang 
[1972]  showed  that  this  assumption  gave  better  accuracy  than  the  assumption  of  plane 
stress  i.e.,  <r33  =  0.  For  a  laminate,  (57)  through  (59)  were  directly  used  to  set  up  the 
constitutive  equations  for  the  first  order  as  well  as  the  higher  order  theories. 

It  is  known  from  elasticity  solutions  [e.g.  Pagano  1969,  1970]  that  the  transverse 
shear  stress  distribution  is  close  to  parabolic  over  the  thickness  of  each  layer.  However, 
the  kinematic  assumptions  (37)  and  (38)  correspond  to  shear  strain  independent  of  x3 
within  the  layer.  If,  in  addition,  the  interlayer  continuity  of  tractions  is  enforced,  the 
shear  stress  distribution  would  have  to  be  constant  over  the  entire  thickness  of  the 
laminate.  In  general,  the  in-plane  stress  components  are  calculated  using  the  constitutive 
equations  (57)  and  the  remaining  components  are  derived  from  the  stress-equilibrium 
equations  for  the  three-dimensional  configuration.  Also,  direct  use  of  the  constitutive 
relations  as  derived  above  will  not  allow  for  shear  and  transverse  deformation  properly 
and  may  result  in  erroneous  plate  stiffness.  It  has  been  customary  to  use  a  stiffness 
factor  to  account  for  this  discrepancy.  As  part  of  the  present  research  effort,  consistent 
theories  to  allow  for  shear  and  transverse  deformation  were  developed.  This  is  briefly 
described  in  Section  V.  Detailed  discussion  is  available  in  items  B.1.2,  B.2.6,  B.2.7,  B.3.3, 


B.4.2  and  B.4.5  of  Appendix  B. 


SECTION  V 

A  VARIATIONAL  APPROACH  TO  CONSTITUTIVE 
THEORY  OF  LAMINATED  PLATES 


5.1  PRELIMINARIES 

Reissner  [1984]  derived  consistent  shear  deformation  constitutive  equations  from  a 
mixed  variational  principle  proposed  on  a  somewhat  ad  hoc  basis.  Herein,  we  start  by 
writing  a  general  variational  theorem  for  elastic  bodies  following  Sandhu  [l 975,  1977] 
and  then  specialize  it  to  the  case  of  a  laminated  plate.  Reissner's  functional  is  seen  to 
arise  as  a  special  case  of  the  general  formulation.  Specialization  to  the  case  of  a 
laminated  plate  using  an  equilibrium  stress  field  yields  the  coupled  constitutive 
relations  between  force  resultants  and  laminar  deformations. 

The  field  equations  of  three-dimensional  linear  elasticity  written  as  a  self-adjoint 
operator  matrix  in  the  complementary  formulation  [Sandhu  1975]  are 


(68) 


Here  8a„  is  the  identity  tensor,  -3-  =  and 

di  ax, 

l  = —( 8  +  §  a  ^ 

1  =  4+6  f  J-) 

2  2  v  “>  $8  nl‘  Qy  ' 


(69) 

(70) 


35 


He  we  have  separated  the  ’’in-plane"  displacements  and  stresses  from  the  transxerse 
and  shear  components.  If  the  bilinear  mapping  is  the  inner  product 


f.gdR 


(71) 


where  f,  g  are  functions  defined  oxer  the  region  R,  dements  of  the  operator  matrix  in 
(68)  satisfy  self-adjointness  in  the  sense  of  (A.2l)  of  Appendix  A  to  this  report.  The 
operators  on  the  diagonal  are  symmetric  and  the  off-diagonal  operators  constitute  adjoint 
pairs  as  defined  in  Appendix  A.  Consistent  Ixrundary  conditions  associated  with  the 
field  equations  are 

u  n  =  u  n  on  S,  (72) 

■  j  1  j  i 

— cr  n  =  — t  on  S,  (73) 

ji  J  i  2 

Here  a  superposed  circumflex  denotes  the  prescribed  value  of  the  quantity  over  the 
portion  of  the  boundary  indicated  in  the  equation,  t,  are  the  cartesian  components  of 
the  traction  vector  and  n,  are  components  of  the  outward  normal  to  the  boundary.  S, 
and  S2  constitute  complementary  subsets  of  the  boundary  of  the  solid.  Besides  the 
boundary  conditions  listed  above,  there  could  additionally  be  discontinuity  conditions  on 
the  variables  at  surfaces  interior  to  R.  These  can  be  stated  as 

(unj  =  g  on  S.  (74) 

—  (a  n )  =  —  h  on  S,  (75) 

ij  i  j  2i 

Here  the  superscripted  prime  denotes  a  jump  in  the  quantity  and  Sls,  Sa  are  contained 
in  the  union  of  surfaces  in  the  interior  of  R.  If  there  is  no  delamination  and  there  are 
no  discontinuities  in  tractions  in  the  physical  problem,  the  right-hand  sides  of  (74) 
and  (75)  vanish,  i.e.,  the  continuity  condition  may  be  regarded  as  a  "homogeneous” 
discontinuity  equation. 
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5.2  A  VARIATIONAL  PRINCIPLE. 


Using  the  definition  (A. 20),  the  governing  functional  for  the  field  equations  of 
linear  elasticity  along  with  the  consistent  boundary  conditions  and  internal 
discontinuity  conditions  is 

fi(u  ,<r  )  =  <u  ,  <r  +  2f  >,  +  «r  ,  —  u  +  C, ,  or.  +  <cr  ,un  — 2un  >c 

r  ij  j  jjtj  j  K  ij  ».J  k. l»j  kl  K  ij  i  j  i  j  Sj 


—  <u  .  <7  n  —  2t  >c  (76) 

1  J1  J  lb. 

It  is  easy  to  show  that  the  Gateaux  differential  of  this  functional  along  arbitrary 
paths  vanishes  if  and  only  if  the  field  equations  along  with  all  the  consistent 
boundary  conditions  and  internal  discontinuity  conditions  are  satisfied.  Noting  that,  for 
sufficiently  smooth  field  variables,  the  following  identity  holds: 


fu  a  dR  =  -  fu  a  dR+  f  ua  n  dS+  |  (ua  n' 

j  i  ji.j  4  ,,J  IJ  J  '  11 1  4  ’  " J 


dS 


(77) 


where  S  is  the  boundary  of  R  and  S,  is  the  set  of  internal  surfaces  imbedded  in  R 
over  which  the  quantity  (u.ovnp  is  discontinuous.  Extended  forms  of  the  governing 
functional  (76)  can  be  written  through  elimination  of  one  or  the  other  of  the  operators 
appearing  on  either  side  of  the  equality  (77).  Elimination  of  derivatives  of  all  the 
components  of  stress  leads  to  a  generalization  of  the  Hellinger  Reissner  functional,  viz., 

Q  ,(u  ,a  )  =  2<u  ,f  >  +  <cr  2u  +  C. ,  cr,  ,>  +  2<cr  ,u  n  -un>  +  2<u  ,t  > 

1  1  ij  i  i  R  ij  i.J  KHJ  kl  R  i.l  i  j  1  J  s  1  1  s 


+  2<cr  ,(u  n  )  —  g  >  +  <u  ,h>  (78) 

■J  1  J  b,J  sjl  1  1  s2l 

Specialization  of  this  functional,  to  the  case  where  the  in-plane  strains  are  derived 
from  stresses,  leads  to  an  extension  of  the  functional  recently  postulated  by  Reissner 
[1984],  Assuming  that  the  boundary  conditions  on  S,  are  identically  satisfied,  the 


functional  reduces  to 


=  £(  2ui  f,  +  O’ja  (  -  2  U3.3  +  ^33}  +  2Cro3(  U  3,q  “  Ua,3  +  KJ  +  ^afi  (  ~2  +  %}) 


dR 


+  2 


f  ut.dS  +  2  f  uhdS  +  2  f  a  jin  u} -gJdS 

O.; 


(79) 


Here  we  have  again  separated  some  of  the  in-plane  quantities  and  introduced  symbols 


e  „  =  C,,  „cr,  _  +  C  ,  .o'  .  +  C  ,  „<r  , 

afi  33a/3  33  y3a/5  >3  y6a/3  yc 


e  =  C,,  ,cr, ,  +  2  C  ,  ,cr  ,+C,  ,<x  . 

o3  33a3  33  y3a»3  >3  >603  yo 


e33  ^3333CT33  +  2  ^>333°"y3  +  ^>f  33°"yfc 


(80) 

(81) 

(82) 


Vanishing  of  the  Gateaux  differential  along  arbitrary  paths  denoted  by  (v,,^)  gives, 
upon  dropping  the  factor  2  throughout, 

A,  ,11,=  f  (t  (  — u  „  +  e  „)  +  t  ,(  —  u  ,  —  u  +  2 e  )  +  t  (  — u  +e  ))dR 

rij,vi'  ^  I  v  a/3  a,/3  ap  a3  a, 3  3,o  a3  33  3,3  33  ' 


-  £ +  °-03(vo.3  +  v3.«)  +  ff33V3.3  “  V,f,)dR 

+  fv,dS  +  f  v.h_dS+  f  TIJ((nJU|)’  —  g,j)dS  +  f  crinj(vi)'dS  =  0 
*2  hi  hi  hi 

Using  the  identity  (77)  and  noting  that  v,  =  0  on  S, 


(83) 


Cl.-  f (v(ov ..  +  {.)  + t  X— u  +e  J  +  t  ,(— u  —  u  +2e  )+t  (— u  +e  ))dR 

(r^v.;  2  I  v  i  j*»j  1  <*/3  <*.p  a3  a3  3'°  »3  33  3,3  33 


J  vi  J  ~  hi)  dS+  J  Ti/(n  jui>  -  8,p  dS 


(84) 


Vanishing  of  the  Gateaux  differential  for  arbitrary  v,,riJ  implies  equilibrium  over  R 
and  S2  ,  as  well  as  the  constitutive  relations  and  the  discontinuity  conditions.  For  the 
special  case  where  the  components  of  stress  cr^  are  restricted  to  satisfying  exactly  the 
relations 
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5.3  EQUILIBRIUM  STRESS  STATES 
53.1  Preliminaries 

Equation  (85)  through  (87)  can  be  satisfied  by  an  infinite  number  of  stress  states. 
Two  different  approaches  were  studied  in  the  course  of  the  present  research.  These 
are: 

1.  Pointwise  or  ’’strong”  satisfaction  of  equilibrium,  i.e.,  choice  of  stress  states  such 

that 
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on 


R 


(90) 


cr  +  f  =0 
j>.j  1 

In  this  case  (85)  would  be  satisfied  regardless  of  the  choice  of  the  displacement  field. 

2.  Global  or  ’’weak”  satisfaction  of  equilibrium,  i.e.,  assuming  the  components  to  be 
polynomials  in  the  coordinate  x3  such  that  for  a  given  choice  of  the  displacement  field, 
the  disequilibrium  force  vector  cr^+fj  is  orthogonal  to  over  R,  t  —  tj  is  orthogonal  to 
Vj  on  S2  ,  and  t_—  h;  is  orthogonal  to  on  S21.  This  will  imply  satisfaction  of 

equilibrium  of  generalized  force  quantities. 

Under  (l)  above,  two  alternatives,  one  based  on  the  first  order  theorv  and  another 
based  on  a  higher  order  theory  were  studied.  For  (2)  one  case  assuming  displacements 
to  correspond  to  the  first  order  theory  but  the  shear  stresses  to  be  cubic  in  x3  was 
examined. 


5.3.2  A  First  Order  Theory 

Assuming  the  in-plane  stress  components  to  be  linear  in  the  x3  coordinate  and 
proceeding  to  integrate  the  equations  of  equilibrium  it  can  be  shown  [items  B.1.2  and 
B.4.2  Appendix  B]  that,  for  vanishing  body  forces, 

<r(i!  =  LQ(k)  +  £,T(k~°  +  £,Ttlc>  (91) 

or  3  3lxa  *2  or  ^3  o 

Here  Tj1'  is  the  value  of  the  stress  component  <ru3  at  the  top  of  the  kth  layer,  and 


f  =±(i-02) 

'  4 

c2  =  -  J-0  +  302 

-I  +  9  +  r’ 


(92) 

(93) 

(94) 


and 


0  = 


(95) 
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5.3.3  A  Higher  Order  Theory 

Assuming  the  in-plane  stress  components  to  be  quadratic  in  the  x3  coordinate,  using 
definitions  (50)  (52),  and  (56)  specialized  to  the  kth  layer,  it  can  be  shown  [items 
B.2.6,  B.2.7  and  B.4.5  of  Appendix  B]  that 

(k)  _  (k)  M(k)  (k) .  .(k)  (k)  D(k)  ( 

“oft  ~  *>.  +  Mo,S  +  1)3  PaS 


where 


7i(.k)  =  -^-(9  -  6O02) 

•  At 


<k)  12 


77(k)  =  ii(-i  +  ne2) 

< 


Using  equilibrium  equations  (3)  for  i  =  1,  2  along  with  specialization  of  (50)  and  (54) 
to  the  kth  layer,  for  no  body  forces, 

°ik3  =  £(,k)Qlk)  +  $(2k)M +  <3  Vk"l)  +  $(k¥k)  (100) 


where 


^(k)=  J_(l-402) 

J I  Ot 


£>= Mo-4  02) 


£<k)  =  --  +  le  +  3  92  -  io  e3 

J  4  2 


(103) 


<(,k)  =  -4---j0  +  302  + IO©3  (104) 

4  2 

Again,  using  equilibrium  equation  (3)  for  i  =  3  along  with  specialization  of  (55)  to 
the  kth  layer 


<7<k)  =  ftk)  N(k)  +  £(k)  T(k  1 1  +  £(k)T(k)  +  /k)  T<k~1)  +  /k)  T<k) 

33  3 1  33  3  ^3  3  ^4  ^5  p,P 
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(106) 


£(k)=  _i0_(i  -  802  +  16  e4) 

‘  16tR 


£(k)  =  _L( -7-240  +  12O02  +  3203  -  24O04) 


±  =  J_(  -  7  +  24  0  +  1 20  02  -  32  03  -  240  04) 

s3  1 A 


£(!°  =  _!*_(  1  +  8  0  -  24  02  -  32  03  +  80  04) 
4  32 


f4k)  =  -  1  +  8  0  +  24  02  -  32  03  -  80  04) 

5  32 


5.3.4  An  Alternative  Approach:  The  Weak  Form  of  Equilibrium. 

In  this  approach,  equilibrium  is  satisfied  in  the  ’’weak”  sense  i.e.,  for  a  given 
choice  of  the  displacement  feld. 


f (or  .  +  f  )u 

4  J"J  1  ' 


dR  =0 


For  a  first  order  theory,  let  the  displacement  field  be  given  by 

ulk) = vlk)(  V +  x3^1k)(  V  (112) 

U3°  =  V<3k^XP  +  X3k)  (H3) 

Let  the  stress  field  be  such  that  cr^,  cr(Qk',  <7jk)  are,  repectively,  linear,  cubic,  and 
quadratic  in  the  coordinate  x(3k>.  Thus 

%S  =  ai+V(3k)  (II4) 

±3  =  S  +  b2  X3  +  C2  X3  }  +  d2  (X3  }  (l  15) 

' (k)  ,  .  (k)  .  1  (k )\2 

Cr33~a3  +  ^>3X3  +  C3  X3  ^  (116) 

Evaluation  of  coefficients  leads  to  the  following  expressions: 


(k)  (k)  ,  /w  +<k)  -Ik)-, 

u  =  v  +0(u  — u  ) 

00  On 
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(118) 


(k)  (k)  ,  a,  ><k)  (kK 

U3  =V3  +d(U3  U3  } 


(k)  (k)M(k)  ,  (k)w(k) 


ap 


_ (k)  _  Jk)  „(k)  ,  yi. k)  ..(k)  ,  Jk)„(k-1)  ,  Jk)„(k) 

^3"^  +*2  Ma3+^3  To  +  $4  T« 

_ (k)  _  >(k)  KT(k)  >(k)_j(k-l)  ,  »(k)_,(k) 

*33  =  *!  N33  +  ^2  T3  +^3  T3 


where 


(119) 

(120) 
(121) 


(k)  1 

n'  'h 

(122) 

1  k )  _  12/> 

V,  — 

lk 

(123) 

/  =  J_(l-402) 

‘  2tk 

(124) 

^(k)=  30  fl(i_4e2) 

(125) 

/k)  =  _I  +  Ae  +  302-ioe3 

4  2 

(126) 

t[k)  =  -  -  -  -  e  +  3  e2  +  io  e3 

4  4  2 

(127) 

£(k)  =  _1_  (  i  -  4  02) 

‘  2tk 

(128) 

l'k)  =  -  -  -  0  +  3  e2 

2  4 

(129) 

£(k)  =  -I+0  +  302 

3  4 

(130) 

Equation  (ill),  using  (117)  through  (130)  yields  the 

following  equations  of  equilibrium 

for  generalized  forces: 

K,(k)  -*-00  (k) 

N„  „  +  cr  ,  —a  ,  =0 

po,/3  o3  a  3 

(131) 

,.(k>  ,.(k>  .  lk  ,  *(k)  .  (kl) 

M„  n  —  Q  +•— -(cr  +cr  )  =  0 

pn.ft  ^  o3  a  3 

(132) 
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(133) 


~(k)  .  +<k)  -<k)  A 

Q  +°\,  =0 


MU  -  .NT'  +  -^(ff™  +  cr.  )  =  0 

a  3,a  33  2  33  33 


(134) 


These  are  the  same  generalized  force  equilibrium  equations  as  were  obtained  by  direct 
evaluation  of  the  integral  in  (47)  for  no  body  forces  and  n  =  0,  1.  Thus,  the  stress 
state  selected  satisfies  the  equilibrium  equations  in  terms  of  generalized  forces  though 
pointwise  equilibrium  is  not  satisfied. 


5.4  CONSTITUTIVE  RELATIONS  FOR  FORCE  RESULTANTS. 

5.4.1  Preliminaries. 

The  variational  equation  (88)  holds  for  every  state  satisfying  (85)  through  (87) 
along  with  uf  identically  equal  to  the  specified  value  on  the  boundary  S,  and  provided 

that  (77)  is  valid.  Substitution  of  the  expressions  (80)  through  (82)  for  e|k)  and  using 
equilibrium  stress  components  defined  in  Section  5.3  would  give  constitutive  relations 
between  the  force  resultants  and  the  gradients  of  the  displacement  vector.  Using 
polynomial  approximations  (36),  the  constitutive  equations  relating  the  displacement 
parameters  with  force  resultants  are  realized.  Herein,  we  summarize  the  results  of  the 
three  schemes.  The  first  is  based  on  the  first  order  displacement  expressions  (37)  and 
(38)  along  with  the  equilibrium  stress  field  based  on  cr^1  linear  in  the  coordinate  Xjk>. 
The  second  is  based  on  u^0  quadratic  and  u,10  linear  in  x^1  with  the  equilibrium 
stress  field  described  by  (96),  (100),  and  (105).  The  thrid  approach  uses  the  weak 
form  of  equilibrium  described  by  (ill)  along  with  displacement  assumptions  (i!2)  and 
(113)  leading  to  the  stress  state  described  by  (119)  through  (121). 


44 


5.4.2  A  First  Order  Shear  Deformation  Theory. 

Assuming  that  the  in-plane  strains  e^1  =  u(u  c)  are  identically  equal  to  i.e.,  the 
in-plane  stresses  are  defined  by  the  displacement  assumptions  and  the  material 
constitutive  laws,  and  noting  that  for  the  first  order  theory  (38)  implies  633  =  ^3  =  0, 
neglecting  e3j,  the  following  specialization  of  (88)  is  realized  for  the  case  where  S(l  is 
identified  as  the  set  of  internal  surfaces  on  which  a  slip  of  magnitude  gu3  occurs: 


/ 


N  2 


v*  r  -oof  w  <k>  ,  ..-wx , 

>  /  r  I  —  u  —  u,  +2e  ,)dx. 

I  q3v  0.3  3.o 


k=l  q 
2 


d  A  =0 


(135) 


Substituting  (37),  (38)  for  u^k)  and  Ujk)  ,  using  (90)  and  cr^k'  in  (80),  and  carrying  out 
the  integration  over  the  thickness,  for  the  equality  to  hold  for  arbitrary  values  of  r^’ 
satisfying  equilibrium. 


,(k)  (Id  _  .  r- 

d>  +w.  =4C  ,  , 

'Or  3,0  >3o3 


(136) 


6  ~(k)  _  _l_/._(k- 1)  -.(kk 

St.  10  >  > 

k 

g(k)  —  n  u (k)  —  n  uJk)  =  C<k3)  ,(-3Q(k)-t.  T(l  l)  +  4t  Tlk)) 

®or3  3  a  a  3  >3a3v  v y  k  y  k  >  J 

+  C(k  *  —  3  Q(k+ 1 5  —  t  T(k’)}  +  4t  T(k>) 

>3o3v  k+J  y  k+ l  y  y 

for  k  =  1,2 .  N-l  (137) 

These  are  2(2  N—l)  equations  in  as  many  unknowns.  Written  in  matrix  form  for  all 
the  layers,  (136)  and  ( 1 37)  constitute  a  set  of  narrow  band  equations  in  which  the 
shear  deformation  of  any  layer  depends  upon  the  shearing  force  in  that  layer  and  the 


shearing  stresses  at  the  interfaces  above  and  below  that  layer.  Equations  (137) 
represent  displacement  discontinuity  equations  for  the  interface.  In  case  there  is  no 
internal  slip  i.e.,  g(0V  =  0  and  the  displacements  are  continuous,  i.e..  u,‘kl  and  Ujk)  vanish, 
(136)  and  (137)  for  all  the  layers  can  be  written  as 


ka  k,lb 

\a 

R 

T 

— 

a 

+ 

^M>b 

\ 

0 

Tb 
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where 


Xa  is  the  set  {Q(uk>,  k  =  1, 2, . . .  N  },  Xb  is  the  set  {o-^’,  k  =  1,2, ...N-l  },  Ra  is  the  set 
{0lk)  +  uw  k  =  l,2,...N}  and  Ta,  Tb  depend  upon  the  specified  shear  forces  at  the  top 
and  the  bottom  surfaces  of  the  laminate.  Elimination  of  Xb  through  static  condensation 
yields 


‘.(kijf  ,(k)j  j  X 

+V3.P} 


j-1 


(139) 


and,  conversely, 


,oo  . 

</>  +v 


3,o 


(140) 


where  j  denotes  the  layer  number.  The  matrices  with  elements  and  \^)J  are  full 

and  symmetric.  For  a  homogeneous  plate,  (139)  gives 

q. = { >  + + <3  >  <>4» 

where  Ea3S3  =  C„3S3  .  For  no  surface  shears,  this  reduces  to  Reissner's  well-known  shear 
deformation  equation. 

For  the  classical  laminated  plate  theory,  <f>  1!°  =  0a  for  all  k.  Then,  defining 


Q„  =  ZQlk). 


Qo  =  £LC(^+u3.P  (,42) 

k  - 1  j-1 

5.4.3  A  Higher  Order  Coupled  Constitutive  Theory. 

Using  the  kinematic  assumptions  of  u0  quadratic  and  u  linear  in  x,  ,  substituting 
for  strains  in  terms  of  the  kinematic  field  variables,  and  writing  e(J  in  terms  of  force 
resultants  using  (96)  to  (110)  ,  the  variational  equality  (88)  yields  the  constitutive 
relations  for  the  force  resultants.  Explicitly,  the  variational  equality  has  the  form 
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(k)  (k)  (k)  ^(k) 

-<11,  ,V2  .1)3  >cy633 
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(k) 

N33 

T.(k- 1 ) 


(k) 

T3 

lrp(k-  i)| 

1  AP 

y(k) 

1  AP 


■  (k)  /  \ 
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|dx3  =  0  (143) 


Carrying  out  the  integrations  over  the  thickness  and  writing 
y  -  2e  =  ©  +v, 

'  a  a3  r&  3.0r 


(k)  -  . (k)  ,  ,(k) 

K  =  20  +0, 

a  a  r  3,a 


(144) 

(145) 


the  equality  holding  for  arbitrary  values  of  Qlk> ,  MlV ,  T*ak) ,  ,  T,3k)  leads 


to  the  following  constitutive  relations  for  the  case  where  identically  equals 


( uanfl) 
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(146) 


(k)  _  _4_  p(k) 

Ka  7»  >3or3 

k 
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4tk  >6  t3  I'6 


(147) 


+  C*k) 

+  U33aS 
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Equations  (146)  through  (151)  are  the  14  constitutive  relations  for  the  14  force 
resultants  for  each  layer.  Equations  (152)  and  ( 1 53)  are  the  continuity  equations.  If 
interlayer  slip  and  delamination  were  not  present,  the  left-hand  sides  of  these  equations 
would  be  zero. 


5.4.4  An  Alternative  Theory. 

Using  the  kinematic  assumptions  (112)  and  (113),  substituting  for  strains  in  terms 
of  the  kinematic  field  variables  and  writing  e^'  in  terms  of  the  force  resultants  using 
(114)  through  (116)  along  with  the  material  constitutive  relations  (80)  through  (82), 
the  variational  equality  yields  the  constitutive  relations  for  the  force  resultants. 
Explicitly,  the  equality  has  the  form: 
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(154) 


Carrying  out  the  integrations  over  the  thickness,  for  the  equality  to  hold  for  arbitrary 
values  of  the  quantities  ’andT^  1  leads  to  the  following 

constitutive  relations: 
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Equations  (155)  through  (159)  are  11  N  constitutive  equations  for  the  mechanical 

variables  N^,  Q(ak),  and  The  equations  (160)  and  (161)  are  the  displacement 

continuity  conditions.  In  case  there  are  no  displacement  discontinuities  in  the  physical 
problem,  the  left-hand  sides  of  these  two  equations  will  vanish. 

5.5  DISCUSSION. 

A  comparison  of  the  three  alternative  approaches  described  above  shows  that  the 
constitutive  relationships  are  identical  wherever  the  field  variables  involved  are  the 
same.  For  example,  (136),  (146)  and  (157)  are  identical  relationships  for  shear 
deformation.  Similarly,  (147)  and  (158)  for  the  bending  stiffness  are  identical.  Other 

relationships  have  diffeient  form  because  different  field  variables  are  used  to  describe 
the  deformation. 

These  theories  explicitly  allow  for  the  laminated  nature  of  the  system  and  do 

away  with  the  need  for  use  of  ad  hoc  stiffness  factors.  The  relations  are  compliance 
type  i.e.,  they  express  deformations  in  terms  of  force  resultants.  Thus  they  can  be 

directly  included  in  a  complementary  type  variational  formulation.  Direct  formulation 
is  possible  by  inverting  the  constitutive  relations.  Indeed,  for  the  first  order  theory, 
(139)  represents  the  stiffness  relations  for  the  shearing  force.  It  should  be  noted  that 

even  though  the  complementary  relationships  are  narrow  band,  the  stiffness  relations 
will  be  full  matrices  i.e.,  the  force  in  any  layer  will  cause  deformations  in  all  the 

layers.  For  the  higher  order  theories,  inversion  of  the  constitutive  relations  would 
obviously  be  more  complicated.  On  the  free  edge,  some  of  the  traction  components  are 
specified.  These  cannot  be  condensed  out  as  could  be  done  for  the  interior  points.  Here 
the  full  form  of  the  constitutive  matrix  has  to  b-.-  retained.  The  computer 

implementation  of  the  application  of  finite  element  methods  to  laminates  described  in 
Section  VIII  allowed  for  this  feature. 
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SECTION  VI 

A  RESTATEMENT  OF  PAGANO'S  THEORY 


6.1  INTRODUCTION 

Solving  three-dimensional  equations  of  linear  elasticity  for  a  homogeneous 
orthotropic  laminated  plate,  Pagano  [1969,  1970]  observed  that  the  in-plane 

displacements  were  layerwise  almost  linear,  and  the  shear  stresses  in  each  layer  were 
nearly  parabolic.  Following  Reissner's  procedure  for  plates,  assuming  the  in-plane  stresses 
to  be  linear  over  each  layer,  Pagano  [1978]  solved  the  equations  of  equilibrium  to 
obtain  shearing  stresses  which  were  quadratic  and  transverse  normal  stress  which  was 
cubic  in  the  transverse  coordinate.  Thus  all  the  six  components  of  stress  were 
considered  in  the  analysis.  The  constitutive  equations  were  derived  using  a  variational 
approach.  Pagano  satisfied  the  continuity  of  displacements  as  well  as  tractions  at  the 
interlayer  surfaces.  The  accuracy  of  results  from  Pagano's  theory  depends  solely  upon 
the  accuracy  of  the  assumption  of  linear  in-plane  stresses  over  each  layeT.  Thus,  the 
procedure  can  be  made  as  accurate  as  desired  by  subdividing  each  lamina  into  a 
number  of  layers.  However,  the  theory  yields  13N  (where  N  is  the  number  of  layers) 
coupled  equations  which  are  difficult  to  solve.  To  overcome  this  difficulty,  Pagano  and 
Soni  [1983]  proposed  a  ’’global-local”  approach  to  the  problem.  In  this  the  laminate 
would  be  modelled  by  a  higher  order  theory,  e.g.,  Whitney's,  to  get  approximate  values 
of  stresses  and  deformations.  This  result  would  then  be  used  as  input  to  a  second  stage 
solution,  based  on  Pagano's  theory,  for  the  local  problem.  However,  the  global  solution 
being  grossly  in  error,  the  final  solution  is  still  generally  unreliable  [e.g.,  Hsiao  and 
Soni  1987].  Numerical  solutions,  using  the  finite  element  or  finite  difference  methods. 


53 


are  possible.  For  application  of  the  finite  element  method  in  a  systematic  fashion  it  is 
desirable  to  write  a  variational  principle  for  the  problem. 

As  part  of  the  present  research  effort,  Pagano’s  theory  was  carefully  examined. 
The  constitutive  relations  for  the  force  resultants  were  obtained  by  direct  integration  of 
strains  and  their  moments  without  resorting  to  any  variational  principles.  The  number 
of  field  variables  was  reduced  by  elimination  of  the  quantities  and  Consistent 

boundary  conditions  were  identified  and  variational  principles  developed  for  the 
problem.  Extended  variational  formulations  admitting  relaxation  of  the  continuity 
requirements  on  some  of  the  variables  were  introduced  along  with  suitable 
specializations  to  make  the  problem  tractable.  One  of  the  specializations  was 
implemented  in  a  finite  element  computer  program.  Details  of  this  investigation  are 
contained  in  items  B.1.8,  B.1.9,  and  B.4.3  of  Appendix  B.  Herein,  we  merely  present  an 
outline  of  the  principal  concepts  and  summarize  the  resulting  equations.  Variational 
formulation  of  Pagano's  theory  as  well  as  the  displacement-based  theories  described  in 
Sections  IV  and  V  is  discussed  in  Section  VTI. 


6.2  PAGANO'S  THEORY. 


6.2.1  Equations  for  Stress  Components. 

Assuming  the  in-plane  stresses  to  be  linear  over  a  layer,  using  the  midsurface  as 
the  origin  for  the  local  transverse  coordinate  x3,  we  have: 


N  .  x,M  „ 
o-  =-S£  +  12  3 


(162) 


Substituting  this  expression  in  the  equations  of  equilibrium,  Pagano  [1978]  evaluated  the 
remaining  components  of  stress  as: 


tr  ,  =  ~~(Q 
0,3  2t  °  2  n 


2  x,  , 
1  -(^-)2 


+  T  +(x,+  -)  — 

°  3  2  t 


(163) 
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o-  =I(T  +T*)-1(T+  -T  )  —  —  —  —  (T+  4-  T  )h.-h. 

33  2  33  33  ^  °».o  o.a  ^  ^  ofo  o,o  2  4 


_  jL(  j  —  T+)  —  —  — *3 
2'  3  3'  t  ,t3 


where 


P  =T  -T 

Of  o  o 


S  =  — (T+  +  T  ) 
a  2  0  0 


P«g;i no's  results  are  identical  to  Reissner's  though  the  form  is  somewhat  different. 

6.2.2  Equilibrium  Equations. 

Substitution  of  the  expressions  for  components  of  stress  in  (55)  followed  by 
integration  leads  to 


2 

N.  =  1(T^  +  T7)+ i-(T+  —  T~  ) 

33  2  ^  3  j  2  or,o  a.o 


Similarly,  defining 


z 

M33  =  /  CT33X3dX3 


We  have,  upon  appropriate  substitutions, 


3  2 

M,,  =  — —  (  T+  +  T~  )+-L-(t!-T:) 
33  120  a-°  ao  10  3  3 


Pagano  combined  ( 1 69)  with  the  equation  of  transverse  force  equilibrium  to  write  the 
equations  of  equilibrium  for  the  plate  as  (48)  ,  (51)  and 

Q  +  M  +  T7  —  T*  —  —  (T+  +  T  )  =  0  (170) 

^0.0  ^2  33  3  3  ^  0,0  0.0 


Some  manipulation  of  (169)  and  (170)  leads  to  the  following  alternative  equation: 
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(171) 


Q«.a  +  — 


+  5(T  -T*)--(T+  +  T  )  =  0 

J  J  2  t»,o 


This  equation  can  be  solved  for  M33  and  the  relationship  used  to  replace  (169)  above. 


6.2.3  Constitutive  Equations  for  Force  Resultants. 

Pagano  [1978]  set  up  a  mixed  variational  principle  to  derive  the  constitutive 
equations  for  the  problem.  This  required  introduction  of  quantities  defined  as  follows: 


f,f',f)  =  f  -  f(x)dx 

J,  »  t2  t3 


where  f  is  any  function  of  x.  In  the  context  of  a  lamina,  components  of  the 
displacement  vector  as  functions  of  the  transverse  coordinate  x3  are  integrated  with 
appropriate  weights  indicated  in  (l  72).  In  our  statement  of  Pagano's  theory,  we  do  not 
use  any  variational  principle  but,  instead,  directly  integrate  strain-stress  relations  (9) 
through  (11).  Using  (l)  to  define  the  strain  components,  and  carrying  out  the 
appropriate  integration,  we  get  the  following  results. 


2 

/  eo0dX3  2ives 


U(a,/3)  t  ^yba0  ^  yl>  +  ^33a/3  ^33  ^ 


X3e<W5dX3 


U,  Ot  =  (  C  c  n  M  ,  +  C,  ^  n  M-  ) 

^2  yoor/3  yo  33  tvp  33 


(174) 


2 

J*e33dx3  £ives 
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(175) 


ik  —  ii,  —  C  r  N  +  C  N 

3  3  o/333  op  3333  33 

_t 

-> 

/  X3e33dX3  g‘VeS 

t 


U3  =  U3+U3“7(C. 


M  +  C  M 

.033  «0  3333  33 


(176) 


f  xje^dx,  gives 


6  u , 


■('..*33-V 


:3333  (T3  +T3 


7(:3333  N33  +  ^(U3"U;) 


/  XJe33dX3  8»VeS 


(177) 


U3'“J(U3+U3  5t(  o/i33Mo^  7(^3333  ^*33  j  (>5  ('3333  *T3  "*  3  ' 


(178) 


/ 

/e„dx3  2ives 


u,  =  -  ('  Q,  —  — (u+  —  u  ) 

3.o  ^  r  3  o  3  ^  p  t  o  o 


(179) 


/  X3eo3dX3  £'VeS 


lu;  --u  =  —  c„,  ,(t^-t:)--(u+  +  u  ) 

4  3,o  2  n  ^  p3cv3  3  3  2  0  0 


(180) 


f  *k  ,dX3  ^'VeS 


u  =  Ac,  (T"+Tj+— c  ()  +  Au  -  Afu"  -  u") 
Vo  l  3i»3  3  3  ^  /’3o3  t  <J  t  °  “ 


(181  ) 


57 


Elimination  of  the  surface  displacements  from  pairs  of  equations  leads  to 


3G,-u  =4rC  „  M  „  +  -2-  C33  M  —  2-  C  ,(T,  —  T  ) 

3  3  a£33  o/3  7^  3333  33  3333  3  3 


6u  =  2  C  fl„N  „  +  i^C„„N„'-C„„(Tt  +  Tj 

3  a/333  ap  5  3333  33  5  3333  3  3 


4.8 

u.  -u,  —  —  u  = - 

3, a  3,o  t  o  15 


C«.3(Tl+TI,-f  C».3<3, 


(182) 

(183) 

(184) 


Eliminating  N33  and  M33  using  (55)  and  (168)  ,  the  constPuti.v  equations  (173)  and 
(174)  become 


lr,v  =1U<  „3  =  —  C  . N  +-C„  Xt!  +  T,  + -(T+  -T  )) 

2  2  y  2  t  /'P0/*  PP  2  33o/3V  3  3  ^  P.P  P./>  y 

•J-r,?  =-u'  =  —  c  „ m  +  J_(il(T3-r:)  +  T+  +t“  ) 

2  2  r  t  (c“'^  ^3  Mr>aP  ap  io  v  t  3  3  p.n  p.pj 


(185) 


186) 


where 


and 


r  =8  -3-  +  S  -9- 

2  pp  qp  Py  Qp 


v  alu 

p  2  p 


Equation  (184)  can  be  written  as 


(187) 


(188) 


v  =  2(u  -G  )  +  2u‘  =  -L(12C  ,  3Q  —  t (T3  +  TT)) 

3  ,P  4  3,p  3,p  t  p  5t  >3p3  >  3  3 


where  we  have  introduced  the  definitions 


(189) 


"T  3  * 

<p  =  — u 
^p  t  p 


(190) 


(191) 


Equations  (185)  ,  (188)  along  with  (182)  and  (183)  are  the  constitutive  equations  of 
Pagano's  theory. 
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6.2.4  Interface  Displacements. 

The  displacements  at  the  surfaces  of  the  plate  can  be  obtained  from  (177)  through 
(181)  .  Thus 

u+  =  — t(|-u,  —  4-  u  —  -^-u)  —  (—  u  —  -i-  u  ) 
a  8  3'a  8  3“  2t  “  4  3a  2  ° 

+  <192) 

u  =t(— u  —  —  u  —  —  u  )  — (— u  —  —  u  ) 
a  8  3,0  8  3t°  2t  °  4  3,0  2  a 

"93> 


»1  -  7(5iJ  -  “■>> + f  ui + + TP  - 7t  xj  - 30  M„) 


uI = -r«a3  -  sp  -  f  u;  -  %i(<!  (6t; + rp  -  n  n33  + »  m33) 


6.2.5  Continuity  Conditions. 

Continuity  of  tractions  and  displacements  across  the  interfaces  between  the  kth  and 
the  (k+l)th  layer  of  a  laminate  simply  implies 

<3k)=<r0  d96) 

lVl)  (197) 


6.3  THE  FIELD  EQUATIONS  IN  OPERATOR  MATRIX  FORM. 


6.3.1  Equilibrium  and  Constitutive  Relationships. 

The  equations  of  equilibrium  (48),  (51)  and  (171)  along  with  the  constitutive 
equations  (182),  (183),  (185)  and  (189)  can  be  written  collectively  in  the  form 

[A](k){u},lt)-i-[B]<kV}'<k)  +  [C]<k){crr<l',  =  0  for  k  =  1,  2 .  N  (198) 
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and 


(199) 


(200) 


(201) 


(202) 


t>0 


{u} 10  = 


6.3.2  Displacement  Continuity  Equations. 

Substitution  of  (192)  through  (195)  into  (196)  gives 

[A j'k,{<rl  tk  1  +  [B]lk){u}<k)  +  <k>  +  [C]<k  n {u}lk+ ' }  +  IA] <k  *'  {<?}  (k* 1 '  =  o 
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and 
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^'3333tk*l 

dydP 
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6.3.3  The  Forcing  Functions. 

The  quantities  {cr}*m  and  (cr}<N>  are  known  for  a  problem  and,  therefore,  can  be 
moved  to  the  right  side  of  the  equalities  (198)  and  (204)  for  k  corresponding  to  the 
top  and  the  bottom  layers.  These  will  appear  as  forcing  functions  for  the  problem. 
Defining  these  as  [Qll),  [qTn",  [P]"1,  [PjN’,  we  have 
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[Qf,)  =  lAf,W0,= 
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In  (210)  through  (213)  the  superposed  hat  denotes  the  specified  value  of  the  quantity. 
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bA  DISCUSSION. 


The  equations  of  Pagano's  theory  stated  above  are  self-adjoint  in  the  sense  of 
(A.21).  To  complete  the  statement  of  the  boundary  value  problem,  we  need  to  identify 
consistent  boundary  operators.  This  feature  is  discussed,  along  with  the  broader  problem 
of  variational  formulation  of  the  governing  equations,  in  Section  VII. 


64 


SECTION  VII 


VARIATIONAL  FORMULATION  OF  THEORIES  OF 
LAMINATED  PLATES 


7.1  INTRODUCTION. 

In  developing  the  constitutive  theory  for  laminated  plates  in  Section  V,  care  was 
exercised  to  assume  the  path  of  the  Gateaux  differential  to  lie  in  the  same  space  as 
the  one  in  which  the  solution  for  stresses  lies.  This  resulted  in  constitutive  relations, 
for  the  force  resultants,  which  are  symmetric  and  a  strain  energy  function  or 
complementary  energy  function  can  be  written  for  the  plate.  Similarly,  the  equations 
of  Pagano's  theory  are  self-adjoint  in  the  sense  of  (A.21).  In  both  the  assumed 
displacement  and  the  assumed  stress  approaches,  the  complete  set  of  field  equations  can 
be  written  in  the  form 

[X]  {y}  =  {z}  (214) 

where  [X]  is  the  matrix  of  field  operators,  {y}  is  the  set  of  field  variables,  and  {z}  is 
the  set  of  forcing  functions.  These  field  equations  include  the  equations  of  equilibrium, 
the  constitutive  equations  as  well  as  the  continuity  conditions.  For  each  of  the  theories 
described  in  Sections  V  and  VI,  the  system  of  coupled  equations  is  self-adjoint.  If 
consistent  boundary  conditions  can  be  identified,  general  variational  formulations  along 
with  appropriate  extensions  and  suitable  specializations  can  be  systematically  generated 
following  the  procedures  outlined  in  Appendix  A  and  detailed  in  references  [Sandhu 
1970,  1971,  1975,  1976].  Detailed  discussion  of  application  to  laminated  plates,  using 
different  theories  as  the  starting  point,  is  given  in  items  B.1.1,  B.1.2,  B.1.8,  B.4.2,  and 
B.4.3  of  Appendix  B.  Herein,  we  merely  present  an  outline  of  the  essential  features. 
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12  THE  FIELD  EQUATIONS. 

In  (214)  the  operator  matrix  [X]  consists  of  sets  of  operators  applicable  to  each 
layer.  If  the  field  variables  are  arranged  as 


{cr}^'0 

{u}(N) 

then  the  matrix  [X]  has  the  form: 


[A/U  [Bf0  0  0  0  0 

[E [sf0  [Cf2*  [Xf2)  0  0 

0  [C]<2)  lAf2)  [Bf2)  0  0 

o  [a/2)  [b 12)  [sP  [Cf3)  [aP 

0  0  0  [C]'3)  [Ap  [bP 

0  0  0  [A]'3’  [Bf31  [^3> 

0  0  [cP 

0  0  [Ap 


.  [cfN,>  [cfN_,)  [xp1’  0 
o  IaT”  IaT0  Irf"-0  0 

0  0  [SP0  isP”  icP 

0  [cfN)  [aP 
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This  operator  matrix  includes  the  equations  of  equilibrium  as  well  the  constitutive 
equations  for  each  layer  and  also  the  continuity  equations  across  interfaces.  Explicitly, 
the  equations  of  equilibrium  and  the  constitutive  equations  have  been  stated  in  (198), 
those  of  displacement  continuity  in  (205),  and  (210)  through  (213)  define  the  forcing 
functions  for  each  layer.  Collectively,  the  forcing  function  {z}  has  the  form: 

— {P}(° 

-<Q}(1) 

0 

0 

0 

0 


{z}=  0  (217) 

0 

0 

0 

0 

0 

— (q}<n1) 

-{P}(N) 

In  (216)  the  operators  [A?*0  for  k  =  1,2 .  N  and  [E]<k)  for  k  =  1,  2  N-l  are 

self-adjoint  with  respect  to  inner  product  and  the  sets  of  operators  {[A]<k,,[A]<k)}; 
{[dk>,[au};  {[BTK) , [Blk1}  constitute  adjoint  pairs.  This  is  sufficient  to  satisfy  the 
self-ad  jointness  condition,  described  in  Appendix  A  to  this  report,  proposed  by  Sandhu 
[1976]  for  linear  coupled  problems. 
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7.3  CONSISTENT  BOUNDARY  CONDITIONS. 


Consistent  sets  of  boundary  conditions  must  be  identified  for  systematic 
development  of  variational  principles  for  a  given  set  of  field  equations.  For  tha  case  of 
laminated  plate  theories  described  herein,  the  form  of  the  boundary  conditions  was 
found  to  be 

[W]{y}  =  {g}  (218) 

where 
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M(,) 

{*}-<» 

<u}(2) 

{u}(3) 

{a}"<3) 


i  i<i) 

<gj(2) 

<g(r}(2) 

<gj(3) 

{gj(3) 


{y}  = 


and  {g}  = 


(220) 


{u}(N  O 
{cr}  (N  0 


kj 

<gj 


(N-l) 


(N  - 1 ) 


(N) 


In  this  set,  the  operator  [Df10  describes  the  boundary  conditions  for  the  field  variables 
(displacements  and  force  resultants)  for  the  kth  layer.  The  other  set  of  equations,  viz., 

[d]fk_,>{o-r(lt_l> + [Efk) {u}(k) + i#Vr(k)Ff +,) {u}(k*,) + = {gjU)  (22D 


denotes  the  displacement  continuity  equations  at  the  boundary.  The  right-hand  side  of 
(218)  lists  the  specified  values  of  the  field  variables  or  their  appropriate  functions  at 
the  boundary.  Detailed  form  of  indvidual  operators  is  given  in  appropriate  items  listed 
in  Appendix  B.  As  an  illustration,  we  list  herein  the  operators  as  they  appear  in 
Pagano's  theory. 
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7.3.1  CONSISTENT  BOUNDARY  OPERATORS  IN  PAGANO'S 
THEORY 

The  operators  and  quantities  appearing  in  the  vectors  of  field  variables  and  forcing 
functions  in  (219)  and  (220),  for  Pagano's  theory,  are: 


/  i  ( k. )  _  3 

lul  “ 


[Dfk)  = 


0  0  0  -7)  0  0 

0  0  0  0  -7).  0 

0  0  0  0  0  -7) 

T)  0  0  0  0  0 

0  rip  0  0  0  0 

0  0  rip  0  0  0 


,k)  ooo  t x'^L-r)  — dk)„7)  o 

[ETk  =  12  k  “^33  y  10  a^33  > 


(224) 


0  0  0 


o  o  o  -Lttc(k’„7i  — clk'  tj  o 
f FT k  —  12  k  “033  '  10  0^33  'y 


1  ,4k) 


0  0  0 


0  0 


_  1  3,4  k)  _g_  _J3_  2-^k) 

[0](k)  =  140  k  3333^>Qy  420  k  3333^y 

0  0 

— _i_t3r(k)  -n  JL  _li_t2r(k> 
[0fk)=  140  k  3333^p  Qy  420  k  3333^y 


1  r.3i-,k)  L  .3  /^k4l)u.  9  11  r .2  ,4k*l)_  ,2,4k) 

1  _  ins  T|'^3333  +  tk4'k'3333^pgy  210  tk41^'3333  tk^333337| 
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and 


{cr}  u'  = 


(k)  ^>3 


Jump  discontinuities  in  the  field  variables  can  be  introduced  as  conditions  on 
internal  surfaces.  Indeed,  as  the  finite  element  bases  generally  have  limited  smoothness 
across  interelement  boundaries,  even  in  cases  where  there  are  no  discontinuities  in  the 
physical  problem,  the  smoothness  condition  of  the  physical  problem  needs  to  be 
introduced  as  a  set  of  ’’homogeneous”  discontinuity  conditions.  Similar  to  the  format 
tor  the  conditions  at  the  external  boundary,  the  discontinuity  conditions  are: 


tu 

8, 

'(k) 

83 

<k) 

85 

(v(Qk))\ 

•<k) 

82 

(k) 

84 

<*2%. 

•tv) 

86 

l-rzpK.’i  ^  + 1  -  '*> 


.  , — .( k )  r  1  /x.(k)y  ,  1  /u(kKi 

+  C  V  l~— (N  J  +  — (M  JJ 

a^33  >12  10  aP 

.  1  ,  3  ,4k- I )  .  3,4k)  x  ,  <k)4 

+ - ( t,  t.C..,,,  Jr)  (<r  .  ) 

k-1  3333  k  3333  'y  p).P 

11  /  2  ,4k-l)  2,4k)  x  ,  (k)x'  ,  1  p/k+l)  ,v(k^l)4 

+  2T^(tk-.C3333  -lkC3333)  75y(<r33  )  +  T2tkMCo033^y(‘\0  } 

,  1  ,4k^ 1)  _  i>V  I  ,4k*  D  ,  <k-ih' 

+  JO  o033  ^Mo0  140tk*l<"3333  T,>  (Crp3,p 

+  —  t.2  C(k:,,T)  (cr,<1k‘,))  =gtk)  on  S(k)  (231) 

420  k  1  3333  'y  33  btr  i 

The  subscripted  i  indicates  an  internal  surface.  The  "homogeneous"  discontinuity 
condition,  i.e.,  vanishing  g,"’1,  represents  the  internal  continuity  constraint  for  the 
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7.4  A  VARIATIONAL  PRINCIPLE. 

}'or  the  set  of  field  equations,  including  the  continuity  equations,  along  with  the 
boundary  conditions  and  the  jump  discontinuity  conditions  given  by  (214),  (218),  (230), 
and  (231),  a  variational  principle  can  be  constructed  following  the  procedure  outlined 
in  Appendix  A.  The  governing  functional  for  the  problem,  following  (A JO),  is: 

n({u},(cr})  =  2<{u}(1),  {P}(1)>s(11  +  2<{u}(N>}(N)>r(n, 

+  £*Iu},k,.[cf  V  <k  '  >r(U+  £<{u}(k,,[A]<k,{u}(k)>R(k) 

k  2  k  I 

+  £<fu}(k\(BlkV}  (k'v, 

k  I 

+  r<{cr)  <k/,[B](k){u}(k)>  <u 

k  =  l  k= 1 

+  £<{cr}'(k),[Cfk+,){u}(k*1)>R(k^)+  £<{tr}‘(k),[Afk){cr}'<k^ ”>R(k) 

k« 1  k=2 

k  =  I 

+  2<l<r}  IN  ".[of  +  2<(<r>  "Mof 1 

+  £<(ur.[D(u(u|,l'-2lgJ>“, 

k  1 

+  [Hf 1 W 11 + '  W ' 1  ’ + 1  W2) + iefV} (2)  -  2{gj 1 

+  £<{cr}  <k) , [efV)  <k  °  +  [Ff k){u}(k)  + 

k-2 

+ [rf- 1 W"  ” + [fff  "Vi""-2{gor',»sii( 

+  <i<ri K  ".[of  "m  ,N  »+i0f  "(uf  ,l+[*f  "Vi"N  ,,+[n,N,w'v 

i'N  I) 

'S(N  I) 
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Gateaux  differentiation  for  various  theories  developed  in  the  course  of  the  present 
research  program  are  contained  in  appropriate  items  listed  in  Appendix  B. 
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7.5  EXTENDED  VARIATIONAL  PRINCIPLES. 

In  the  operator  matrix  [X]  in  (216)  the  sets  {[Bjki  ,[5]“°};  {IC]“‘\[C]'k,k  and 
{[A]10 .[A]00}  constitute  pairs  of  adjoint  operator  matrices.  Thus 

<{u}(k> . [cfk){oT(k_,)>  =  <{cr}  <k  11 , 


(k  1)  „(k)  fk^ik),  1  -.(kk 

>3  'y  o/333  12  a&  10  S°° 


(233) 


~(k  1)  -r-tk)  (k) 

<cr  ,A  cr  > 


(k)  .(kK-(k-l) 


,  —  «J  ,  A  cr 
<  * 


D(k) 


,  _  (k  I >  1  3„(k)  (k)  (k)  1  3,.(kl  _  (k  1) 

+  <CT  ,  T)  ,  - t,C,,,,CT  ,  — <cr  T)  , - 1,  ,  ,0"  ,  >  ,,  , 

>3  'y  J4()  k  3333  p3,p  s<k>  y3  'y  }  4()  k  3333  n3.p  s'k) 

tk)  13  2„(k)  ~  (k  1 )  .  ~  (k  1  >  13  ,2,Jk>  ~(k> 

+  <a%,  - t.C.,,,,0-  ,  n  >...  +  «r,,  ,  — - 1  .  cr  ,  ti  > 

33  420  k  3333  >3  'y  kk)  33  k  777  »  \.3  >■ 


42()  k  3333  y3  ''y  s' 


(k> 


(234) 


and 


<{u}(k) .  [Bf  k){oT(k)>  =  <{cr}“(k\[B](k){u}<k)> 

,  _-(k)  ^(k)  r  1  M^k)  *k  *.(kk 

+  <tr>3  ^y  ,Co^33^  1Q  12  Na^s(k) 


(235) 


Also,  the  operator  matrices  [Ef0  and  [Af*5  are  self-adjoint.  This  implies  that  some  of 
the  operators  in  these  matrices  can  be  replaced  by  their  adjoints  with  appropriate 
changes  in  the  boundary  terms.  Thus,  it  is  possible  to  eliminate  the  operators  requiring 
differentiation  of  the  generalized  force  quantities.  This  extends  the  class  of  admissible 
functions  to  include  force  resultants  which  may  not  be  differentiable.  Specifically, 
elimination  of  operators  A^’,  A^1,  and  A^’  from  the  operator  matrix  [AT*’  can  be 
accomplished  through  the  following  identities: 


—  (k)  »(k).,(k)  v*k)  A(k)_(k)  »,\iw  ~\« 

<v  ,A,.N  n>  ,.,  =  <N  ,A.,v  >,.,  +  <N  ar)  ,v. 

y  14  op  p(k)  pp  41  y  j^vk)  op  ’&  p 


.(k)  _(k) 


c(k) 


,(k)  » (k)» *<k)  x.(k)  A<k),<k)  .  ».(k)  ,(k) 

r  25  np  u<k>  pp  52^  jj(k)  ap  'a  j(k) 


~(k)  A(k),.(k)  ,.(k)  _(k) 

<v  ,  A,,Q  >  ..  ,  =  <Q  7) 

3  36 'o  R'k*  'o  ro  3  ^'k) 


(236) 

(237) 

(238) 


74 


Using  (236)  through  (238)  to  eliminate  terms  containing  derivatives  of  force  resultants 
leads  to  functionals  with  extended  domain  of  definition.  Some  of  these  functionals  can 
be  specialized  to  yield  formulations  in  which  the  number  of  field  variables  is  reduced 
by  requiring  some  of  the  field  equations  to  be  satisfied  identically.  One  such 
specialization  is  stated  below. 


7.6  A  SPECIAL  VARIATIONAL  PRINCIPLE. 

Assuming  the  variational  formulation  to  be  extended  through  elimination  of 
derivatives  of  the  force  resultants,  and  further  assuming  that  the  physical  problem  has 
no  displacement  and  force  discontinuities,  i.e.,  (gu) k’  vanish,  and  that  the  displacement 
boundary  conditions  are  satisfied,  the  governing  functional  reduces  to  the  one  used  by 
Chyou  [items  B.1.8  and  B.4.3,  Appendix  B]  to  set  up  a  finite  element  solution  scheme 
for  the  problem.  This  functional,  in  explicit  form,  is: 


r\  f  x  — 1*1  10)  «  -Ai)  *  _(o)  , 

ni(u,o-)  =  2<v>  ,ayi>  +<$y  ,t1or>3>  +2<v3  ,a 


-(1)  _<o) 


rvi'  y  1  yj  k 

-  — (N)  (1st)  ,  -AS)  .  (N)  ,  -j  -<N)  __(N> 
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k=  2 
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ofi  10  <*£33  k  >3  33  l^(k)  oP  jQ  afl33  >3,y  ^  33  j 

k  1  k 


„(k)  2  ~(k)  -<k  1)  , 

-"Qp  *  5C>3p3°'>3  >)} 


N 
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N  1  t 

.  — <k)  (k)  .  -r<k)  lk  „  (k)  ,  -(k)  _  <k)  « 
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(239) 


Because  the  mechanical  variables  M^’,  and  Q^k)  are  completely  defined  by  the 
constitutive  relations  (185),  (186)  and  (189),  the  above  specialization  has  only  the  field 
•  ariables,  v^’,  <J*uk),  v^’  and  crf3k>.  The  number  of  field  variables,  therefore,  is  8N, 
where  N  is  the  number  of  layers.  This  is  less  than  the  13N  in  Pagano’s  original 
theory.  Other  specializations,  allowing  for  elimination  of  only  some  of  the  terms 
containing  derivatives  of  the  force  resultants,  are  possible.  Alternatively,  terms 


76 


containing  derivatives  of  the  kinematic  'ariables  can  be  eliminated  using  the 
self-adjointness  of  the  field  operator  matrix.  The  particular  specialization  described  by 
(239)  was  used,  in  the  present  research  program,  to  develop  finite  element  solution 
schemes  for  Pagano's  theory  as  well  as  the  assumed-displacement  layerwise  theories  of 
laminated  composite  plates.  Detailed  discussion  of  the  finite  element  implementation  is 
given  in  relevant  items  listed  in  Appendix  B. 
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SECTION  VIII 


FINITE  ELEMENT  STUDIES 

8.1  INTRODUCTION. 

The  theoretical  development  described  in  the  previous  chapters  was  used  to  develop 
computational  models  for  stress  and  deformation  analysis  of  laminated  composites.  The 
work  consisted  of  three  groups  of  studies  based  on  the  finite  element  method,  viz.. 

1.  Prliminary  studies  directed  towards  accurate  determination  of  stresses  from  dis 
placement  solutions  determined  by  finite  element  methods. 

2.  Elastic  analysis  of  free-edge  delamination  specimens. 

3.  Analysis  of  laminated  composites  including  free-edge  delamination  specimens 
using  laminated  plate  theories. 

In  this  section  we  briefly  describe  each  of  the  three  components. 

8.2  PRELIMINARY  STUDIES. 

The  finite  element  method  for  analysis  of  linear  elastic  solids,  based  on 
minimization  of  potential  energy  using  displacement  interpolation  functions  which  do 
not  have  continuous  derivatives  across  element  boundaries,  yields  discontinuous 
approximation  to  the  stress  field.  Stresses  across  interfaces  as  well  as  at  the  surface  of 
the  solid  are,  therefore,  difficult  to  determine.  To  ensure  stress  continuity  across 
element  boundaries,  Hinton  and  Campbell  [1974]  proposed  a  smoothing  function 
approach.  A  continuous  stress  pattern  over  the  entire  domain  was  obtained.  Oden  and 
Brauchli  [1971]  derived  improved  stresses  based  on  the  theory  of  conjugate 
approximations.  However,  the  determination  of  the  conjugate  functions  involved  the 
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inversion  of  a  rather  large  matrix.  An  approximate  implementation  of  the  concept, 
using  consistent  nodal  averages,  was  proposed  by  (Men  and  Reddy.  Zienkiewicz  et  al. 
[1974]  proposed  use  of  a  least  squares  procedure  with  a  mixed  formulation  to  arrive  at 
a  continuous  stress  approximation.  Mirza  and  Olson  [1980]  studied  the  performance  of 
the  mixed  finite  element  method.  Chang  [1981  ]  and  Jenq  [1982]  tried  several  different 
interpolation  schemes  and  found  that  these  procedures  were  expensive  and  not  always 
sufficiently  accurate.  Moazzami  [1984]  extended  Cook's  work  [1982]  on  Loubignac's 
method  to  plane  elasticity.  This  method  gave  good  approximation  for  stresses  but  was 
quite  expensive  because  of  increased  bandwidth.  Cook  and  Huang  [1986]  proposed  that 
instead  of  fitting  the  average  nodal  stress  to  the  equilibrium  equation  over  each 
element  by  using  Loubignac  iteration,  an  approximate  nodal  strain  that  made  use  of  a 
finite  difference  scheme  be  used.  In  the  present  research  effort,  a  scheme  for 
calculation  of  nodal  point  stresses  using  the  method  of  undetermined  coefficients  applied 
to  the  displacement  solution  from  a  reliable  finite  element  analysis  was  considered. 
Interpolation  between  the  nodal  points  would  define  the  continuous  field.  The  method 
was  applied  to  several  plane  elasticity  problems  for  which  the  exact  solutions  are 
known.  A  four-point  isoparametric  quadrilateral  element  modified  to  eliminate  the 
spurious  shear  modes  was  used  to  get  the  displacement  field.  Two  approximation 
schemes,  one  based  on  making  the  formula  exact  for  quadratic  polynomial  interpolation 
for  displacement  and  the  other  for  cubic  polynomials,  were  used.  Averaging  procedures 
based  on  different  selections  of  the  nodal  points  in  any  group  used  for  the 
determination  of  the  formula  were  tried.  It  was  found  that  the  procedure  could  yield 
good  estimates  of  stresses  even  for  relatively  coarse  meshes. 

Another  study  considered  the  use  of  higher  order  elements  with  continuous 
gradients  across  element  boundaries.  Such  elements  were  originally  introduced  by  Clough 
and  Felippa  [  1 968]  for  analysis  of  plate  bending.  Tocher  and  llartz  [1967]  extended 
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their  application  to  finite  element  analysis  of  plane  elasticity  and  showed  the 
effectiveness  of  ’’conforming”  elements.  In  the  present  research.  Tocher  and  Hartz'  work 
was  revisited.  Clough  and  Felippa's  triangular  elements  were  combined  into 

quadrilaterals.  Figure  3  shows  Felippa's  LCCT-12  cubic  displacement  triangular  element. 
The  quantities  wy  and  w2  denote  derivatives  of  the  field  variable  w  with  respect  to 
the  spatial  coordinates  y  and  z,  respectively.  Assuming  normal  gradients  to  vary 
linearly  along  an  edge  the  midside  nodal  point  on  that  edge  can  be  eliminated.  This 
results  in  elements  designated  as  LCCT-11  or  LCCT-9  depending  upon  whether  one  or 
all  three  of  the  midside  nodal  points  are  eliminated  in  this  manner.  LCCT-9  is  the 

element  used  by  Tocher  and  Hartz  [1967].  Figure  4  shows  quadrilaterals  built  up  from 
four  LCCT-12  elements  and  its  variants,  the  LCCT-11  and  the  LCCT-9.  These 

quadrilaterals  are  referred  to  as  Q-23,  Q-19,  and  Q-15,  respectively.  For  homogeneous 
materials,  these  elements  yielded  stress  distribution  close  to  the  exact  even  for  very 

coarse  meshes.  Figure  5  and  Figure  6  show  finite  element  models  of  two  example 
problems  solved.  Figure  7  through  Figure  10  show  the  results  of  the  finite  element 
analyses  compared  to  the  Euler-Bernoulli  beam  theory  in  one  case  and  a  Ritz  solution 
obtained  by  Hiremath  [1985]  in  the  other.  Results  from  the  conventional  four-point 
isoparametric  element,  designated  Q- 4,  are  also  included.  In  the  present  research  the 
Q-23  element  was  developed  for  analysis  of  laminated  composites.  However,  in  that 
context,  it  has  to  be  kept  in  mind  that  the  tractions,  and  not  strains,  are  continuous 
across  the  inter-laminar  boundaries. 
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Figure  4: 


(a)  LCCT-12  (b)LCCT-ll  (c)  LCCT-9 


(a)  Q-23  (b)  Q-19  (c)Q-15 


QUADRILATHRAI.  HFFMFN’TS:  Q-23,  Q-19,  AND  Q-15. 
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FINITE  ELEMENT  MODELS  ()1  SIMPLE:  BE.AM. 
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Figure  7:  SIMPLE  BEAM:  Y-DISPLACEMENT  AT  POINT  A. 
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Figure  8:  SIMPLE  BEAM:  LONGITUDINAL  STRESS  AT  POINT  A. 
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Figure  9:  PLATE  UNDER  PARABOLIC  LOAD:  X-STRESS  AT  POINT  A. 
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Figure  10: 


PLATE  UNDER  PARABOLIC  LOAD:  X-STRESS  ALONG  X=0. 


8.3  ANALYSIS  OF  FREE-EDGE  DELAMINATION  SPECIMENS. 

8.3.1  Introduction. 

Success  of  the  continuous  strain  finite  element  model  in  solving  problems  of  plane 
elasticity  motivated  its  extension  to  the  problem  of  free-edge  delamination.  Specialization 
of  the  three-dimensional  elasticity  problem  to  the  case  of  free-edge  delamination 
specimens  results  in  three  components  of  displacement  depending  upon  only  two  spatial 
variables.  For  this  pseudo-two-dimensional  problem,  the  continuous  strain  element  can  be 
used  with  only  minor  changes  to  yield  a  fully  three-dimensional  stress  state.  However, 
in  the  actual  specimen,  the  tractions  are  continuous  at  interlaminar  surfaces  but  not  so 
the  strains.  The  continuous-strain  representation  would  result  in  discontinuous  tractions. 
In  the  present  research,  this  observation  led  to  the  development  of  a  continuous-traction 
finite  element  representation.  This  procedure  involves  cubic  interpolation  of  the 
displacement  components  over  each  element  and  can  be  quite  expensive  to  use  for  large 
problems  involving  a  large  number  of  layers.  The  Air  Force  had  previously  used  an 
axisymmetric  model  to  represent  a  free-edge  delamination  specimen.  This  is  because  an 
annulus  of  large  internal  radius,  subjected  to  a  uniform  internal  radial  pressure  of 
small  magnitude,  will  stretch  essentially  uniformly  in  the  tangential  direction.  Thus,  a 
small  segment  along  the  circumferential  direction  will  be  a  practically  straight  coupon 
of  uniform  width  under  uniform  longitudinal  strain.  In  the  present  research  effort, 
this  procedure  was  also  revisited.  The  method  was  extended  to  include  a  complete 
three-dimensional  state  of  stress  by  permitting  circumferential  axisymmetric  disp'  cement 
in  addition  to  the  radial  and  transverse  displacements. 

In  the  following  sections,  we  describe  the  finite  element  implementation  of  the 
continuous  strain  and  continuous  traction  finite  element  analyses  as  well  as  the 
axisymmetric  modelling  of  free-edge  delamination  specimens. 
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8.3.2  Continuous  Strain  Finite  Element  Modelling  of  Free-Edge 
Delamination  Specimens. 

For  plate  bending  analysis,  Clough  and  Felippa  [1968]  used  the  values  of  w,  the 
transverse  displacement  of  the  plate,  and  its  derivatives  w,  and  wy  as  the  generalized 
coordinates  in  the  finite  element  representation.  For  the  free-edge  delamination 
specimen,  the  displacement  field  involves  three  components  of  displacement  dependent 
on  two  spatial  coordinates.  Thus,  nodal  values  of  u,  uy,  u2;  v,  vy,  v2  would  be  required 
in  addition  to  w,  wy,  and  w?.  Derivation  of  the  interpolation  functions  is  given  in 

items  B.1.5,  and  B.3.1  of  Appendix  B.  Figure  11  shows  the  nodal  point  degrees  of 

freedom  for  the  element  LCCT-12  for  the  pseudo-two-dimensional  elasticity  problem. 

Elements  of  quadrilateral  shape  can  be  set  up  as  assemblages  of  triangular  elements  as 
discussed  in  Section  8.2.  The  element  Q-23  has  23  degrees  of  freedom  for  each 
component  of  displacement.  However,  the  interior  nodal  points  can  be  eliminated  by 
static  condensation  to  yield  an  element  with  48  degrees  of  freedom  represented  by  the 
values  of  the  three  components  of  displacement  along  with  their  derivatives  at  the 
four  corners  and  the  normal  derivatives  of  the  displacement  components  at  the  four 
midside  nodes.  Application  of  the  procedure  to  cross-ply  and  angle-ply  free-edge 
delamination  specimens  showed  [item  B.3.1,  Appendix  B]  that  the  element  could  provide 
good  estimates  of  stress-distribution  over  the  specimen  except  for  the  component  cryz. 
This  was  because  the  element  could  not  satisfy  the  traction-free  condition  at  the  free 
edges.  This  was  the  motivation  for  the  development  of  the  continuous  traction 

procedure  described  in  the  next  section. 
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8.33  Continuous  Traction  Finite  Element  Modelling  of  Free-Edge 
Delamination  Specimens. 

High  interlaminar  stresses  near  the  free  edge  are  primarily  responsible  for  the  onset 
of  delamination  in  free-edge  delamination  specimens.  In  order  to  use  stress-based 
criteria  for  onset  and  propagation  of  damage,  it  is  necessary  that  a  reliable  estimate  of 
interlaminar  stresses  be  available  for  a  given  situation.  The  continuous-strain  element 
discussed  in  the  previous  section  is  useful  for  homogenous  materials  but  cannot  give 
correct  results  for  layered  materials  because  in  these  systems,  it  is  the  interlayer 
tractions  and  not  the  strains  that  are  continuous.  For  the  continuous-strain  element, 
the  components  of  displacement  and  their  derivatives  at  the  nodal  points  were  selected 
as  the  degrees  of  freedom.  Assuming  that  the  interlayer  surfaces  will  lie  along 
element  boundaries,  it  is  reasonable  to  still  use  the  continuous-strain  representation 
within  the  element.  However,  for  connection  with  the  adjoining  elements,  the  quantities 
associated  with  the  nodal  points  should  be  the  traction  components.  If  the  y-axis  is 
parallel  to  the  interfaces  in  a  laminate,  the  displacement  being  continuous  across  the 
interface,  derivatives  with  respect  to  y  would  be  continuous.  Thus  the  derivatives  of  u, 
v,  and  w  with  respect  to  the  z-coordinate  need  to  be  replaced  by  the  components  of 
traction  on  the  interface.  For  nodal  points  located  on  the  interfaces  corresponding  to  z 
=  constant,  the  components  of  traction  equal  the  stress  components  crxz,  ayi,  and  cr^.  For 
midside  nodal  points  and  nodal  points  on  other  inter-element  boundary  orientations,  the 
components  will  equal  the  stresses  orrx,  crnn,  and  crns.  This  requires  a  transformation 
relationship  between  the  degrees  of  freedom  to  be  used  for  ’’global”  assembly  and  those 
to  be  used  for  "local”  or  intra-element  interpolation. 
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83.3.1  Transformation  for  the  Corner  Nodal  Points. 

The  stress-strain  relationship  for  an  orthotropic  material  with  reference  to  the 
cartesian  reference  frame  described  by  x,y,z  axes  is: 
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Here  CtJ  are  components  of  the  compliance  tensor,  written  in  the  reduced  notation,  for 
monoclinic  materials  with  symmetry  with  respect  to  the  x-y  plane.  These  components 
are  related  to  the  material  compliance  coefficients  CtJ  through  the  following 

relationships: 
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where 
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and  m  =  cos  0,  and  n  =  sin  0.  Here  0  is  the  angle  between  the  ply  orientation  (the 
material  1-axis)  and  the  system  x-axis.  Ekk,  G,jt  vtJ,  are  moduli  of  elasticity  in 
extension,  the  shear  moduli,  and  Poisson's  ratios  in  material  coordinates. 

Replacing  the  strains  by  appropriate  expressions  in  terms  of  gradients  of 
displacement  components,  allowing  for  the  special  case  of  derivatives  along  the 
longitudinal  axis  (the  x-axis)  being  zero,  (240)  gives: 


>1  Pi.  D. 2  S 

=2  °21  D22  a2 


where 


{e  }  =  v  If  (  =  w  +  v 
'  r  y  ’  lc2*  y  z 


<T 

z 

,  {o-2}  = 

<T 

yz 

cr 

xz 

c  r 

1  '“'12  S6 

c  c 

2  S2  S6 

,  [D 

c  r 

6  Sb  Sb 

C3b  0  0 


94 


C,3 

c23 

£3* 

£33 

0 

0 

[D;,]  = 

0 

0 

0 

,  ID  J  = 

0 

C44 

C4s 

0 

0 

0 

0 

C4S 

Css 

Elimination  of  {<r,}  between  the  two  equations  in  (243)  gives 
{€*}  =  [C*]{<T2} 

where 
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Solving  (251 )  for  u2,  v2,  and  w2  in  terms  of  the  remaining  quantities  gives  the 
following  transformation  from  the  global  degrees  of  freedom,  which  include  components 
of  traction  at  the  nodal  points,  to  the  displacement  degrees  of  freedom  which  are  used 
to  describe  the  variation  of  displacement  within  the  element: 
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83.3.2  Transformation  for  Midside  Nodes. 

The  transformation  from  stress  components  normal  to  the  element  boundary  at  the 
midside  points  and  the  normal  displacement  gradients  used  for  interpolation  within  the 
element  is: 
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1  i 

y\,  A  are  the  coordinates  of  the  corner  nodal  points  of  the  triangular  finite  element  and 
Stj  are  elements  of  the  material  stiffness  matrix  [S]  which  is  inverse  of  the  compliance 
matrix  [C]  with  elements  CtJ. 


8.3.3.3  Traction-free  Boundary  Conditions. 

The  traction-free  boundary  conditions  for  one  quadrant  of  the  free-edge 
delamination  specimen  are: 

o-Jy,  H)  =  cryz(y,  H)  =  crjy ,  H)  =  0  (267) 

at  the  top  surface  along  with 

or  (B,  z)  =  cr  (B,  z)  =  a  (B,  z)  =  0  (268) 

yy  *-y 

at  the  lateral  free  edge.  Here  2B  and  2H  are  the  total  width  and  the  total  thickness 
of  the  specimen.  The  Q-23  element  has  the  tractions  at  the  top  surface  included  in  the 
list  of  field  variables  for  the  finite  element  model.  Therefore,  these  conditions  can  be 
explicitly  satisfied  simply  by  specifying  the  values  at  the  nodal  points  located  on  the 
top  surface.  Similarly,  for  the  midside  nodal  points  on  the  lateral  free  edge,  the 
tractions  are  field  variables  for  the  finite  element  model  and  their  values  can  be 
specified  directly.  At  the  corner  points  of  the  elements,  the  component  cryx  can  be 
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directly  specified.  However,  the  components  cr  and  tr  have  to  be  specified  as  linear 
constraints  on  the  problem.  Expressing  these  in  terms  of  the  kinematic  variables,  we 
have: 


°’yy  =  0==  S12e0  +  S23Wz  +  S22Vy  +  S26Uy 

=  0  =  S.6e0  +  §36Wz  +  §26Vy  +  §66Uy 
Solving  for  v  and  u  in  terms  of  the  remaining  quantities: 


V 

y 

*„ 

*  1  2 

eo 

u 

y 

*2, 

*22 

w 

z 

where 


^.  =  ^S26S,6-S665I2) 


*12  a  ^26^36  ^66^23^ 


*21  a^S26S12  S22S16^ 


(269) 

(270) 


(271) 


(272) 


*22  a^26S23  ^22S36^ 

and 

ft  =  S22S66-SL  (273) 

Using  (251)  to  write  wz  in  terms  of  cr^,  (271)  can  be  incorporated  into  the 
transformation  (258). 


8.3.4  An  Axisymmetric  Model. 

8.3.4. 1  Introduction. 

A  computer  program  previously  used  by  the  Air  Force  to  model  free-edge 
delamination  specimens  was  applicable  to  an  axisymmetric  configuration.  The  model 
assumed  the  tangential  shearing  strains  to  be  zero  i.e.,  the  deformation  is  purely  radial 
and/or  axial.  It  was  based  on  the  use  of  a  bilinear  Lagrange  interpolation  or  optional 
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use  of  a  reduced  integration  procedure  proposed  by  Singh  [1977]  and  Sandhu  and  Singh 
[1978].  In  the  present  research,  the  program  was  enhanced  to  include  ’’twist”  in  the 
analysis.  The  enhanced  version  has  the  capability,  which  the  earlier  version  lacked,  to 
model  the  complete  stress  state  in  angle-ply  specimens.  Herein,  we  give  an  outline  of 
the  theoretical  development.  Details  are  available  in  items  B.1.6  and  B.1.7  of  Appendix 
B. 


8.3.4.2  The  Axisymmetric  Model. 

The  equations  of  equilibrium,  in  cylindrical  coordinates,  for  linear  elastic 
deformation,  are: 
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The  strain-displacement  relationships  in  cylindrical  coordinates  are: 
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(274) 


(275) 


where,  r,  z,  and  9  are,  respectively,  the  radial,  the  axial  and  the  tangential  coordinates. 
ur,uz,and  ue  are  the  displacements  corresponding  to  the  three  coordinate  axes.  For 
axisymmetry,  all  quantities  are  independent  of  6,  and  the  equations  (274)  and  (275) 
reduce  to 


100 


8.3.5  Computer  Implementation. 

The  schemes  described  in  Sections  1.3.3  and  1.3.4  were  implemented  in  finite 
element  computer  programs.  The  codes  were  verified  by  application  to  Pagano's  problem 
[l  978]  and  then  applied  to  analyze  stresses  and  deformations  in  selected  22-layer 
delamination  specimens. 
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8.3.5.1  Continuous  Traction  Finite  Element  Analysis. 

The  finite  element  discretization  of  a  quadrant  of  the  cross-section  of  a  free-edge 
delamination  specimen  is  shown  in  Figure  12.  Each  layer  was  modelled  by  four 
sublayers.  In  the  transverse  direction,  the  elements  were  of  width  B/10  in  the  middle 
part  and  B/50  for  the  ten  elements  nearest  the  free  edge.  The  cross-ply  specimen  was 
[0/90^  and  the  angle-ply  specimen  was  [  ±45],.  Comparison  of  results  from  the  con¬ 
tinuous  traction  model  with  Pagano's  theoretical  solution  is  given  in  Figure  13  through 
Figure  20.  The  agreement  is  excellent. 

The  program  was  also  used  to  solve  a  22-layer  problem.  Several  different  layups, 
for  which  test  data  were  available,  were  examined.  Details  are  given  in  items  B.1.5  of 
Appendix  B.  Here,  we  show  the  results  for  [( ±49.8)s/90l.  specimen.  Figure  21  shows 
the  finite  element  model  used  for  one  quadrant  of  the  specimen.  Figure  22  shows  the 
distribution  of  cru  along  the  mid-plane  of  the  specimen  in  the  vicinity  of  the  free 
edge.  The  ’’constant  strain  element”  denotes  results  of  a  pseudo-two-dimensional  finite 
element  analysis,  ignoring  trxi  and  crX)  and  using  Q1  elements,  reported  by  Sandhu  and 
Sendeckyj  [1987].  Through-the-thickness  distribution  of  cru  and  crxi  at  the  free  edge  is 
shown  in  Figure  23.  The  plot  shows  considerable  irregularity  in  the  distribution  of 
<tzj  at  the  free  edge.  Figure  24  shows  the  influence  of  distance  from  the  free  edge  on 
through-the-thickness  distribution  of  The  irregularity  in  stress  distribution  is 

confined  to  the  free  edge  and  the  distribution  is  much  smoother  even  a  short  distance 

away  (-^-=0.995)  from  the  free  edge.  The  distribution  of  cr  at  the  free  edge  wras 
B 

carefully  examined  by  refining  the  mesh  in  the  vicinity.  Finite  element  models,  based 
upon  refinement  along  each  of  the  y  and  the  z  axes  and  along  both  the  axes,  were 
used.  Along  the  y-axis,  the  refinement  consi'ted  of  subdividing  each  edge  element  into 
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FINITE  ELEMENT  MODEL  OF  FREE-EDGE  DELAMINATION  SPECIMENS. 
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Figure  13: 


DISTRIBUTION  OF  X-STRESS  ALONG  CENTER  OF  TOP  LAYER  (45 
DEGREES)  IN  A  [45/-4S]  SYMMETRIC  LAMINATE. 
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Figure  14: 


DISTRIBUTION  OF  XY-STRESS  ALONG  CENTER  OF  THE  TOP  LAYER 
(45  DEGREES)  IN  A  [4S/-45]  SYMMETRIC  LAMINATE. 
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238  (14x17)  F.  E.  Mesh  -  Continuous  Traction  Analysis 


Figure  21:  FINITE  ELEMENT  MODEL  OF  22-LA YF.R  DELAMINATION  SPECIMEN. 
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Figure  24:  THROUGH  THE  THICKNESS  DISTRIBUTION  OF  ZSTRESS  NEAR  THE 

FREE  EDGE  OF  THE  22-LAYER  FEDS. 
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four.  Along  the  z-axis,  refinements  only  for  the  90-degree  ply  near  the  midplane  of 
the  specimen  as  well  as  for  all  the  layers  were  carried  out.  Selected  results  are 
plotted  in  Figure  25  and  Figure  26.  Refinement  along  the  z-axis  over  the  90°  layer 
improved  the  solution  over  that  layer  without  much  change  in  the  stress  patterns 
elsewhere.  Refinement  over  the  y-axis,  on  the  other  hand,  improved  the  results  over 
the  layers  other  than  the  90'  layer  with  little  effect  on  the  stresses  in  that  layer.  To 
improve  stress  solution  all  along  the  free  edge,  it  was  necessary  to  refine  the  mesh 
along  both  the  axes.  Figure  27  and  Figure  28  show'  plots  of  the  through-the-thickness 

variation  of  crIZ  at  -^-=0.99.  Refinement  of  the  mesh  has  little  effect  on  the  stress 

distribution  at  this  location.  Thus,  the  influence  of  mesh  refinement  is  entirely 
confined  to  the  immediate  vicinity  of  the  free  edge. 

Figure  29  shows  the  variation  of  <Jyi  along  the  center  of  the  90-degree  layer  for  a 
distance  of  eight  ply  thicknesses  from  the  free  edge. 

The  computer  program  is  quite  efficient.  For  the  154-element  mesh  with  513  nodal 
points,  the  total  number  of  algebraic  equations  was  2311  with  a  band-width  of  186. 
The  execution  time  for  this  case  was  121  seconds  on  an  IBM  3081  mainframe 
computer. 
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Figure  25: 
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HROUGH  THE  THICKNESS  Z-STRESS  DISTRIBUTION  AT  Y/B  =  .99 
DR  THE  22-LA YF.R  FEDS.  (REFINEMENT  ALONG  Y-AXIS  ONLY). 
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8.3.5.2  The  Axisymmetric  Model. 

For  the  axisymmetric  model,  the  modified  four-point  isoparametric  element  was 
used.  The  modification,  introduced  by  Doherty  et  al.  [1969]  and  Zienkiewicz  et  al. 
[1972]  consists  of  evaluating  the  contribution  of  shearing  strains  to  the  energy  using 
one-point  Gauss  quadrature.  This  amounts  to  the  assumption  of  constant  shearing  strain 
in  the  element.  Use  of  this  ’’selective”  reduced  integration  procedure  has  been  found  to 
give  better  estimates  of  stress  in  thin-walled  structural  systems.  Comparison  of  this 
procedure  with  Sandhu  and  Singh's  [1978]  reduced  integration  scheme  showed  that  the 
results  from  the  two  analyses  were  comparable  but  the  ’’selective"  scheme  was  the 
more  economical  in  computational  effort.  To  ensure  uniformity  of  extensional  strain 
over  the  delamination  specimen,  different  values  of  the  internal  radius  of  the  annulus 
were  tried.  It  was  found  that  using  an  internal  radius  exceeding  lO4  times  the 
thickness  of  the  specimen  would  give  essentially  uniform  tangential  strain  in  the 
annulus  corresponding  to  uniform  extension  in  the  specimen. 

Figure  30  shows  the  correspondence  between  the  free-edge  delamination  specimen 
and  a  segment  from  a  large-radius  ring  under  radial  pressure.  A  576-element  mesh 
shown  in  Figure  31  was  used  to  model  the  delamination  specimens  analyzed  by  Pagano 
[1978].  It  should  be  noted  that,  for  the  model  used,  the  complete  thickness  of  the 
laminate  has  to  be  modelled  instead  of  the  half-thickness  used  in  the  continuous 
traction  element  analysis  described  in  the  previous  section.  Comparison  of  results  from 
the  axisymmetric  model  and  the  continuous  traction  element  analyses  showed  that  the 
numerical  results  were  in  close  agreement  while  the  computational  effort  was  much 
smaller  for  the  axisymmetric  model  which  uses  a  lower  order  element.  Of  course,  for 
comparable  accuracy,  a  finer  mesh  would,  in  general,  be  required  for  the  simpler 
model.  As  an  illustration.  Figure  32  shows  through-the-thickness  distribution  of  cr^  at 
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576  (32x18)  Finite  Element  Model  For  Axisymmetric  Analysis. 


Figure  31:  FINITE  ELEMENT  MODEL  OF  THE  AXISYMMETRIC  PROBLEM. 
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Free  Edge 


*  CONTINUOUS  TRACTION  MODEL 

*  32*18 

-  SPLINE  FIT  OF  DATA  INDICATED  BY  * 

WITH  QUADRATIC  FUNCTION  NEAR  INTERFACE 


THROUGH-THE-THICKNESS  Z-STRESS  DISTRIBUTION  AT  Y/B 
FOR  THE  [0/90]  SYMMETRIC  LAMINATE. 


SURFRCE 


—  =0.99  for  a  cross-ply  specimen  over  half  the  thickness,  using  8  equal  subdivisions 


over  the  thickness  of  each  layer  and  18  elements  over  half  the  width  as  shown  in 

Figure  31.  Similar  plots  for  the  stress  component  cryi  at  ^-  =  0.99  are  shown  in 

Figure  33.  Figure  34  through  Figure  36  show  the  through-the-thickness  variation  in 

stress  components  crxz,  azl,  and  cryz,  respectively,  for  the  four-layer  [±45],  specimen  at 

the  plane  -=-=0.99  The  results  shown  were  obtained  using  1152  elements  for  the 
B 

axisymmetric  model  (complete  cross-section)  and  166  elements  for  the  continuous 
traction  model. 

The  procedure  was  applied  to  a  22-layer  delamination  specimen  for  which  test 
data  were  available.  Several  mesh  refinements  were  used.  Figure  37  show's  the  finite 
element  mesh  used  to  get  the  plots  -  Figure  38  and  Figure  39  -  showing 

through-the-thickness  distribution  of  and  cry2  at  ^  =  0.995  for  a  [(±25^/90^ 

delamination  specimen  for  the  axisymmetric  model  as  well  as  the  continuous  traction 
model  discussed  earlier. 
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THROUGH-THE-THICKNESS  YZ-STRESS  DISTRIBUTION  AT  Y/B 
FOR  THE  [0/90]  SYMMETRIC  LAMINATE. 


SURFACE 


+  CONTINUOUS  TRACTION  MODEL 
*  32x18 

-  SPLINE  FIT  OF  DRTfi  INDICATED  BY  * 

WITH  QUADRATIC  FUNCTION  NEAR  INTERFACE 


0J40  o'.  80  1 20  l‘.  60  2.00 


DISTANCE  FROM  CENTER  LINE  Z 


THROUGH-THE-THICKNESS  DISTRIBUTION  OF  Z- STRESS  AT  Y/B 
FOR  THE  [4S/-45]  SYMMETRIC  LAMINATE. 


SURFACE 


XZ  STRESS  X  10E-0S  (PSI) 

.14  -0.08  -0.02  0.04  0.10  0.16 


*  CONTINUOUS  TRACT  I  ON  MODEL 

*  32* 1 8 

-  SPLINE  FIT  OF  ORTR  INDICATED  BT  * 
WITH  QUADRATIC  FUNCTION  NEAR  INTERFACE 


Figure  36:  THROUGH -THE-THICKNESS  DISTRIBUTION  OF  YZ-STRESS  AT  Y/B 

.99  FOR  THE  [45/-4S]  SYMMETRIC  LAMINATE. 
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DISTANCE  FROM  CENTER  LINE  Z/H 


8.4  ANALYSES  BASED  ON  LAMINATED  PLATE  THEORIES. 


8.4.1  Introduction. 

A  more  general  approach  to  the  problem  would  be  applicable  not  only  to  free-edge 
delamination  specimens  but  also  to  plates  of  general  configuration.  For  the  present 
research  effort,  this  involved  implementation  of  appropriate  theories  of  bending  and 
stretching  of  laminated  composite  plates.  To  start,  the  existing  layerwise  theory  of 
laminated  plates  was  implemented  in  a  finite  element  computer  program.  Limitations  of 
this  theory  were  immediately  noticed,  and  subsequent  research  effort  was  directed 
towards  development  of  comprehensive  theoretical  basis  for  the  constitutive  and 
variational  theory  of  laminated  composite  plates  and  its  implementation  in  finite 
element  computer  models.  The  theory  has  been  described  in  Sections  IV  through  VII.  In 
this  section  we  summarize  the  finite  element  implementation  and  list  some  results. 

8.4.2  Discrete  Laminate  Theory. 

Mawenya  and  Davis  [1974]  proposed  a  general  finite  element  formulation  using 
quadratic  isoparametric  multilayer  plate  elements.  However,  they  did  not  give  details  of 
their  theoretical  and  numerical  formulation.  In  the  present  research  program,  the 
existing  theory  of  laminated  plates  was  restated  in  a  self-adjoint  form  and  implemented 
in  a  finite  element  computer  program.  The  program  is  capable  of  using  the  bilinear 
isoparametric  interpolation  as  well  as  biquadratic  Lagrangian  or  8-node  serendipity 
interpolation.  Details  of  this  investigation  are  given  in  item  B.1.1  of  Appendix  B. 
Figure  40  shows  the  three  element  types  available  in  the  program.  Application  to 
sandwich  as  well  as  homogeneous  plates  showed  very  good  performance  in  modelling 
deflections.  However,  results  for  <ra,  crXI,  and  cryi  were  quite  inaccurate  near  the  free 
edge.  It  was  found  that  refinement  of  mesh  oveT  the  thickness  resulted  in 
improvement  in  stress  predictions  while  refinement  in  the  y-direction  did  not  have  a 
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1.  Element  Q4:  (a)Rectangular  parent  (b) Isoparametric  counterpart 


Figure  40:  FINITE  ELEMENT  INTERPOLATION  SCHEMES  USED  FOR  DISCRETE 

LAMINATE  THEORY. 
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significant  effect.  This  was  specially  noticeable  near  the  free  edge.  Figure  41  shows 
the  sequence  of  refinement  of  the  model  along  the  z-axis  i.e.,  over  the  thickness  of  the 
specimen.  Figure  42  and  Figure  43  show  influence  of  refinement  upon  crXI  and  c ryi 
along  the  interface  of  a  [±45],  specimen. 

8.4.3  A  First  Order  Theory. 

The  first  order  theory  for  constitutive  relations  for  force  resultants  in  laminated 
plates  described  in  Section  IV  w'as  written  in  a  self-adjoint  form,  appropriate 
variational  formulations  were  developed  and  a  specialization  was  implemented  in  finite 
element  computer  program.  Items  B.1.2,  11.1.3  B.2.3,  B.2.4,  B.2.6,  B.3.3  and  B.4.2  give 
details  of  the  investigation.  The  finite  element  implementation  used  Hughes'  [1978] 
Heterosis  element  shown  in  Figure  44.  The  traction-free  boundary  condition  could  not 
be  explicitly  satisfied  because  the  interlaminar  tractions  had  been  condensed  out  of  the 
constitutive  equations.  This  condition  could  be  introduced  as  specified  tractions  in  the 
constitutive  theory.  This  would  make  the  constitutive  relations  different  at  the  surface 
and  in  the  interior  necessitating  the  assumption  of  a  pattern  of  variation  from  the  free 
edge  towards  the  inside  of  the  specimen.  Numerical  results  from  the  analyses  showed 
that  the  influence  of  the  shear  coupling  indicated  by  the  theoretical  investigation  was 
not  significant.  This  was  probably  because  the  theory  is  first  order,  based  on  the 
assumption  of  e33  =  0,  and  ignores  the  coupling  between  the  other  force  resultants.  The 
numerical  results  practically  coincided  with  those  from  the  discrete  laminate  theory 
which  did  not  allow  for  the  constitutive  coupling.  It  was  noticed  that  refinement  of 
each  layer  into  sublayers  resulted  in  noticeable  improvement  in  the  solution.  Figure  45 
and  Figure  46  show  the  results  for  subdivision  of  each  of  the  two  layers  in  one-half 
of  a  [±45],  free-edge  delamination  specimen  into  three  and  five  sublayers.  This 
remarkable  improvement  with  refinement  was  the  motivation  for  the  search  for  a 
higher  order  theory. 
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(a)  Two  Sub-layers  (N=2) 


(b)  Six  Sub-layers  (N=6) 


Figure  41:  SEQULNCE  OF  REFINEMENTS  OVER  THE  THICKNESS  OF  ONE-HALF 

OF  THE  [4S/-45]  SYMMETRIC  SPECIMEN. 
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DISTRIBUTION  OF  XZ-STRESS  AT  INTERFACE 
REFINEMENT  OVER  THE  THICKNESS. 
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Figure  46:  DISTRIBUTION  OF  YZ-STRHSS  AT  INTERFACE  OF  THE  [45/-4S] 

SPECIMEN  -  INFLUENCE  OF  REFINEMENT  OVER  THE  THICKNESS. 
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8.4.4  A  Higher  Order  Theory. 

Investigation  of  the  existing  discrete  laminate  theory  and  the  first  order  theory 
showed  that  refinement  over  the  thickness  of  the  laminate  resulted  in  improvement  of 
accuracy  of  the  predicted  stresses.  This  was  the  motivation  for  the  development  of 
higher  order  theories  described  in  Section  VI.  A  specialization  of  the  higher  order 
theory  to  the  free-edge  delamination  specimens  was  implemented  in  a  finite  element 
computer  program  using  cubic  Hermite  polynomial  interpolation  for  all  the  field 
variables.  The  solution  showed  rapid  convergence  with  refinement  of  the  mesh.  Kven 
very  coarse  mesh  (six  elements  over  the  entire  width  of  the  element)  gave  reasonable 
results  for  four-layer  cross-ply  and  [±45l  angle-ply  specimens. 

8.4.5  Pagano's  Theory. 

Pagano's  theory  as  restated  in  Section  VI  was  implemented  in  a  finite  element 
computer  program  using  Heterosis  elements.  Application  to  Pagano's  problems  [1978]  of 
four-layer  cross-ply  and  angle-ply  specimens  showed  excellent  agreement.  The  procedure 
was  then  applied  to  a  22-layer  free-edge  delamination  specimen  [( ±  25.5)s/90\.  For  a 
rather  coarse  64-element  model  shown  in  Figure  47,  the  results  for  cr^  along  the 
midplane  and  along  the  free  edge  are  plotted  in  Figure  48  and  Figure  49.  The  results 
agree  very  well  with  those  from  the  axisymmteric  model  described  earlier. 
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Free  Surface 


256  (4x8x8)  F.  K.  Mesh  -  Pagano’s  Theory 


Figure  47: 


FINITE  ELEMENT  MODEL  OF  THE  22-LAYER  DELAMINATION 
SPECIMEN. 
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o 

o 

o' 


o  Pagano's  Theory  (FKM) 
+  Continuous  Traction  Model 
(no  local  refinement) 


Figure  49:  DISTRIBUTION  OF  Z-STRESS  ALONG  THE  FREE  EDGE  OF  THE 

22-LAYER  SPECIMEN.. 
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SECTION  IX 


DISCUSSION 

For  analysis  of  stresses  and  deformations  in  a  free-edge  delamination  specimen,  the 
continuous  traction  cubic  interpolation  finite  element,  with  appropriate  refinement  of 
the  mesh  in  the  vicinity  of  the  free  edge  and  interfaces,  will  yield  reliable  predictions. 
The  representation  by  a  segment  of  an  annulus  does  not  give  continuous  tractions  but 
is  much  more  economical  in  computational  effort  required  to  obtain  the  solution.  Even 
the  pseudo-two-dimensional  representation  ignoring  axy  and  axz  gives  reasonable  accuracy 
except  near  the  free  edge. 

Existing  theories  of  laminated  plates,  including  the  discrete  laminate  theory,  cannot 
allow  for  the  traction-free  edge  conditions.  Thus,  they  are  not  useful  for  delamination 
prediction  even  though  they  are  adequate  for  representation  of  gross  behavior  like 
deflection  under  load  or  frequency  of  vibration.  Coupled  constitutive  theory,  allow-ing 
for  the  shear  coupling  between  layers,  is  mechanically  consistent  and  independent  of 
the  need  lor  ad  hoc  correction  factors.  However,  it  too  cannot  properly  model  the 
stresses  in  the  vicinity  of  the  free  edge  and  the  interfaces  between  plies.  The  higher 
order  theory,  developed  during  this  reasearch  program,  based  on  quadratic  variation  of 
c ru3  over  the  thickness  of  each  layer  or  sublayer  along  with  use  of  consistent 
constitutive  relationships,  gave  reasonable  accuracy  even  for  coarse  meshes  (six  elements 
over  the  entire  width  of  the  delamination  specimen).  The  computer  program  has  been 
written  for  a  specialization  of  the  theory  to  free-edge  delamination  specimens.  To  solve 
the  problem  of  an  arbitrary  plate  configuration  involving  free  edges,  further  work  on 
development  of  a  comprehensive  computer  program  would  be  necessary.  Finite  element 


implementation  of  Pagano's  theory  shows  that  this  theory,  based  on  the  sole  assumption 
of  linear  variation  of  the  in-plane  stresses  over  each  layer,  is  capable  of  yielding 
accurate  estimates  of  stress  and  deformation  throughout  the  laminate.  As  the  laminae 
can  be  subdivided  into  as  many  sublayers  as  desired,  it  should  be  possible  to  set  up 
sequences  of  solutions  converging  to  the  correct  solution.  However,  this  procedure  is 
quite  expensive  in  computational  effort  as  well  as  storage  requirements.  An  alternative 
solution  procedures  for  this  formulation,  based  on  reduction  of  the  problem  of 
determining  roots  of  the  characteristic  polynomial  to  a  matrix  eigenvalue  problem, 
appears  to  hold  good  promise  of  combining  accuracy  with  economy  of  effort  and 
extending  the  application  of  Pagano's  theory  to  multi-layer  laminates. 
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Appendix  A 

VARIATIONAL  FORMULATION 


Often,  obtaining  an  approximate  solution  to  coupled  boundary  value  problem  relies 
on  appropriate  variational  formulations.  Following  Sandhu's  [1970,  1971,  1975,  1977] 
extension  of  Mikhlin's  [1965]  basic  variational  theorem  to  coupled  linear  boundary 
value  problems  including  nonhomogenous  boundary  conditions,  we  present  a  summary 
of  the  basic  concepts  for  setting  up  the  variational  formulation  applicable  to  the 
problem  of  laminated  plates. 


A.l  PRELIMINARIES 


A.  1.1  Boundary  Value  Problem 

Consider  the  boundary  value  problem 

A  u  =  f  on  R  (A.l) 

C  u  =  g  on  QR  (A.2) 

where  0R  is  the  boundary  of  the  open  connected  region  R  in  an  euclidean  space.  K  is 
the  closure  of  R.  A  and  C  are  the  linear  bounded  operators.  Let  VR  and  VdR  be 

linear  vector  spaces  defined  on  the  regions  indicated  by  the  subscripts,  and  be 

dense  subsets  in  VR  and  VdR,  respectively.  Then  the  differential  A  and  C  can  be 
regarded  as  the  transformations 


A:Wb-Vk 


C:W6^VaK 


(A. 3) 
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A.  1.2  Bilinear  Mapping 

Let  V  and  S  be  linear  vector  spaces.  A  bilinear  mapping  B:VxV-*S  assigns  to 

each  ordered  pair  of  vectors  u ,  v  €  V  an  element  in  S.  Furthermore,  bilinearity  is  sat¬ 
isfied  for  u, ,  u2 ,  v, ,  v3 ,  u ,  v  €  V,  if 

BfaUj  +  u2,v)  =  alKu,  ,v)+  B(u2,v)  (A.4) 

B(u  .oVj  +  v2)  =  aB(u  ,  v()  +  B(u ,  (A.5) 

where  a  is  scalar.  For  convenience,  we  shall  use  the  notation 

BR  (  u  ,  v  )  =  <  u ,  v  >R  ( A.6) 

To  set  up  a  variational  formulation,  symmetric,  nondegenerate  bilinear  mappings  are 

used,  i.e., 

<u,v>R  =  <v,u>R  (A. 7) 

and 

<u,v>R  =  0  for  all  v  if  and  only  if  u  =  0  (A.8) 

A.  1.3  Self-Adjoint  Operator 

An  operator  A*  on  V  is  said  to  be  the  adjoint  of  A  with  respect  to  symmetric 
bilinear  mapping  B^VxV-^S,  where  S  is  a  linear  vector  space,  if 

<  u  ,  Av  >  R  =  <  v  ,  A  u  >  R  +  F>6R  (  v , u  )  ( A.9) 

for  all  u  and  v  €  V  and  where  DdR(u,v)  represents  quantities  associated  with  bound¬ 
ary  QR  of  R.  If  A=  A*,  then  A  is  said  to  be  self-adjoint.  If  A  is  a  self-adpint  opera¬ 

tor,  then  DdR(v,u)  is  antisymmetric,  i.e, 

DeR  (  v  ,  u  )  =  -  DeR  ( u  ,  v  )  (A. 10) 

Furthermore,  A  is  said  to  be  symmetric  with  respect  to  the  bilinear  mapping  if, 

<  u , Av  >R  =  <v,Au>r  (A.l  l) 

The  boundary  operator  C  is  said  to  be  consistent  with  the  self-adjoint  operator  A  if 
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(A. 12) 


DaR(v’u)=  <u’Cv>aR-<v,Cu  >aR 


A.  1.4  Gateaux  Differential 

If  fi:V-*S,  where  V  is  such  that  if  u,  Q6V,  u+  Xu€V  for  scalar  X,  the  Gateaux 
differential  of  a  function  Cl  along  the  path  u  is  defined  by 
ft(u  +  Xu)  —  n(u) 


8  ft(u)  =  lim- 
u  k-0  A 

where  u  is  referred  to  as  the  path. 


(A. 13) 


A.2  THEOREM 


For  the  field  equations  (A.l)  we  def  ine 
ft(u)  =  <u  ,  Au>r  —  2<u  ,  f  >R 

The  Gateaux  differential  of  Cl  is: 

8  n(u)  -jlm  <  u  Xu  ,  A(u  +  Xu)>  —  2<u  +  Xu,f>  —  <u,Au>  +  2<u,f> 
u  U  X 


(A. 14) 


(A. 15) 


=  <  u  ,  Au>  +  <  u  ,  Au >  —  2  <  u  ,f  > 

=  2  <u , Au -f> 

The  Gateaux  differential  vanishes  at  the  solution  u  =  u0  where  Au0  —  f  =  0.  Conversely, 
if  8uft(u)  vanishes  for  all  u,  nondegeneracy  of  <  ,  >  implies  Au0  — f  =  0.  If  the 
range  of  the  bilinear  mapping  is  the  real  line,  vanishing  of  the  function  Cl  would 
imply  its  minimum,  maximum,  or  stationary  value,  depending  upon  the  operator  A 
being  positive,  negative  or  semi-definite. 
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A.3  LINEAR  COUPLED  PROBLEMS 


The  above  discussion  lor  a  single-valued  function  u  can  be  extended  to  the  case  of 
several  variables.  If  there  are  n  variables,  V  is  defined  as  the  direct  sum 

V  =  V .  +  V ,  + . +  V  (A. 16) 

1  2  n 

and  an  element  u6V  is  an  n-tuple  (u,,u2,  —  ,un)  with  u,  €  V,  for  i=l,  2,  n.  A 
bilinear  mapping  on  V  is  defined  as 

<u,v>  =  <u1,v,>Rl+  <u,,v2>r2+ . +  <un,vn>Rn  (A.17) 

where  <  ,  >  is  defined  for  components  ui(v,  of  (uj.lvj  respectively. 

If  the  field  equations  and  the  boundary  conditions  for  a  linear  coupled  boundary 
value  problem  are: 


ZA  u  =  f  on  R 
'j  j  < 


(A. 18) 


j-i 


52  C.j  uj  "  S,  on  0R  i=L2 . 41 

j-i 

the  governing  functional  based  on  Eqs.  (A. 18,  A. 19)  is 


(A. 19) 


n(u)  =  Z  <u,’ZA,Juj_2f.  >R  +  Z<u.’ZcuuJ~2g.> 


SR 


(A.20) 


i*l  j- I  i-l  j-1 

The  set  of  operators  A  is  said  to  be  self-adjoint  with  respect  to  the  bilinear  mapping, 
if 


<  v  ,  7  A  u  >B  =  T*  <  u  ,  A  v  >_  +  D.d(u  ,  v  ) 
1  L*  1J  j  R  Zw  j  ji  i  R  JR  j  I 

r  *  r  i 


(A.21) 


where  DgpCuj.v,)  represents  quantities  associated  with  boundary  QR  of  R.  The  boundary 
operator  C0  are  said  to  be  consistent  with  the  set  of  field  operators  A,j  if 


d.  (u  ,  v )  =  y  <  u .  c  v  >  —  <  v ,  y  C  u  > 

t»R  j  i  £-*  j  j'  '  SR  i  i.i  j 


SR 


(A.22) 


j  i 


162 


Appendix  B 

PUBLICATIONS  AND  PRESENTATIONS 


B.1  RESEARCH  REPORTS 

1.  M.  Moazzami,  R.S.  Sandhu,  and  W.E.  Wolfe,  ”A  Finite  Element  Procedure  for 
Analysis  of  Laminated  Composite  Plates”,  Wright  Laboratory  Technical  Report 
WL-TR-91-3023,  Flight  Dynamics  Directorate,  Wright  Laboratory,  Wright- 
Patterson  Air  Force  Base,  Ohio,  1991. 

2.  R.S.  Sandhu,  W.E.  Wolfe,  SJ.  Hong,  and  H.S.  Chohan,  ”A  Consistent  Shear 

Deformable  Theory  for  Laminated  Plates,”  Wright  Laboratory  Technical  Report 
WL-TR-9 1-3021,  Flight  Dynamics  Directorate,  Wright-Patterson  Air  Force  Base, 
Ohio,  1991. 

3.  R.S.  Sandhu,  W.E.  Wolfe,  SJ.  Hong,  and  H.S.  Chohan,  ”A  Computer  Program 

for  Analysis  of  Laminated  Plates  Based  on  a  Consistent  Shear  Deformable 
Theory,”  Wright  Laboratory  Technical  Report  WL-TR-91-3018,  Flight  Dynamics 
Directorate,  Wright-Patterson  Air  Force  Base,  Ohio,  1991. 

4.  R.S.  Sandhu,  W.E.  Wolfe,  C.C.  Chang,  and  H.R.  Chu,  "Computer  Program 

CTQ23  for  Finite  Element  Analysis  of  Free-Edge  Delamination  Specimens”, 
Wright  Laboratory  Technical  Report  WL-TR-91-3017,  Flight  Dynamics  Director¬ 
ate,  Wright-Patterson  Air  Force  Base,  Ohio,  1991. 

5.  R.S.  Sandhu,  W.E.  Wolfe,  R.L.  Sierakowski,  C.C.  Chang,  and  H.R.  Chu,  "Finite 
Element  Analysis  of  Free-Edge  Delamination  in  Laminated  Composite  Speci¬ 
mens,”  Wright  Laboratory  Technical  Report  WL-TR-91-3022,  Flight  Dynamics 
Directorate,  Wright-Patterson  Air  Force  Base,  Ohio,  1991. 


163 


6.  R.S.  Sandhu,  R.L.  Sierakowski,  W.E.  Wolfe,  and  R.A.  Dandan,  ’’Finite  Element 
Analysis  of  Laminated  Composite  Axisymmetri c  Solids,  Vol.  1  -  Theory  and 
Applications,”  Wright  Laboratory  Technical  Report  WL-TR-91-3020,  Vol.  1, 
Flight  Dynamics  Directorate,  Wright-Patterson  Air  Force  Base,  Ohio,  1991. 

7.  R.S.  Sandhu,  W.E.  Wolfe,  and  R.A.  Dandan,  ”A  Computer  Program  for  Finite 
Element  Analysis  of  Laminated  Composite  Axisymmetric  Solids,”  Wright  Labo¬ 
ratory  Technical  Report  WL-TR-91-3020,  Vol.  2,  Flight  Dynamics  Directorate, 
Wright-Patterson  Air  Force  Base,  Ohio,  1991. 

8.  R.S.  Sandhu,  W.E.  Wolfe,  and  11.11.  Chvou,  ’’Variational  Formulation  and  Finite 
Element  Implementation  of  Pagano's  Theory  of  Laminated  Plates,”  Wright  Lab¬ 
oratory  Technical  Report  WL-TR-91-3016,  Flight  Dynamics  Directorate,  Wright- 
Patterson  Air  Force  Base,  Ohio,  1991. 

9.  R.S.  Sandhu,  W.E.  Wolfe,  and  H.H.  Chyou,  ”A  Finite  Element  Computer  Pro¬ 
gram  for  Analysis  of  Laminated  Plates  by  Pagano's  Theory,”  Wright  Laborato¬ 
ry  Technical  Report  WL-TR-91-3019,  Flight  Dynamics  Directorate,  Wright- 
Patterson  Air  Force  Base,  Ohio,  1991. 

B.2  CONFERENCE  PROCEEDINGS. 

1.  C.C.  Chang  and  R.S.  Sandhu,  ’’Finite  Element  Analysis  of  Free  FJdge  Effect  in 
Composite  Laminates,”  11th  Canadian  Congress  for  Applied  Mechanics,  Edmon¬ 
ton,  Alberta,  June  1987. 

2.  C.C.  Chang,  R.S.  Sandhu,  R.L.  Sierakowski,  and  W.E.  Wolfe,  ’’Continuous  Strain 
Finite  Element  Analysis  of  Free  Edge  Delamination  Specimens,”  ASME  Pressue 
Vessel  and  Piping  Conference,  San  Diego,  California,  June  1987. 

3.  S.J.  Hong  and  R.S.  Sandhu,  ”A  Consistent  Shear  Deformation  Theory  for  Lami¬ 
nated  Composite  Plates,”  13th  Mechanics  of  Composites  Review,  Florida,  Novem¬ 
ber  1988. 


164 


4.  S.J.  Hong,  R.S.  Sandhu,  and  H.S.  Chohan,  ’’Consistent  Shear  Constitutive  Rela¬ 
tions  for  Laminated  Plates”,  12th  Canadian  Congree  for  Applied  Mechanics, 
Ottawa,  Ontario,  June  1989. 

5.  R.S.  Sandhu,  ’’Some  Recent  Developments  in  Mechanics  of  Laminated  Plates,” 
International  Conference  on  Mechanics,  Physics,  and  Structure  of  Materials;  A 
Celebration  of  Aristotle's  23  Centuries,  Thessalonika,  Greece,  August  1990. 

6.  R.S.  Sandhu  and  M.  Moazzami,  ”A  Variational  Approach  for  Consistent  Devel¬ 
opment  of  Constitutive  Relations  for  -Laminated  Plates”,  15th  Mechanics  of 
Composites  Review,  Dayton,  October  24-  25,  1990. 

7.  R.S.  Sandhu,  and  M.  Moazzami,  ’’Constitutive  Relations  for  l-'orce  Resultants  in 
Laminated  Plates,”  Third  International  Conference  on  Constitutive  l^aws  for 
Engineering  Materials:  Theory  and  Applications,  Tucson,  Arizona,  January  1991. 

B.3  REFEREED  JOURNAL  ARTICLES. 

1.  Chern-Chi  Chang,  Ranbir  S.  Sandhu,  Robert  L  Sierakowski,  and  William  E 
Wolfe,  ’’Continuous  Strain  Finite  Element  Analysis  of  Free-Edge  Effect  in 
Laminated  Composite  Specimens,”  American  Society  for  Testing  and  Materials, 
Journal  of  Composite  Technology  and  Research,  Vol.  10,  No.  2,  54-64,  Summer 
1988. 

2.  C.C.  Chang,  R.S.  Sandhu,  R.L.  Sierakowski,  and  W.E.  W’olfe,  ’’Continuous  Stress 
Finite  Element  Analysis  of  A  Free-Edge  Delamination  Specimen,”  Computers 
and  Structures,  Vol.  29,  No.  5,  783-793,  1988. 

3.  Soon  Jo  Hong,  Ranbir  S.  Sandhu,  and  Harpal  S.  Chohan,  ’’Consistent  Shear 
Constitutive  Relations  for  a  Laminated  Plate”,  under  review  for  possible  pub¬ 
lication  in  ASTM  Journal  of  Composite  Technology  and  Research. 


165 


B.4  DISSERTATIONS  AND  THESES. 

1.  Chern-Chi  Chang,  ’’Finite  Flement  Analysis  of  Laminated  Composite  Free-Edge 
Delamination  Specimens”,  Ph.D.  Dissertation,  The  Ohio  State  University,  Colum¬ 
bus,  Ohio,  1987, 

2.  Soon  Jo  Hong,  ”A  Consistent  Shear  Deformable  Theory  for  the  Vibration  of 
Laminated  Plates”,  Ph.D.  Dissertation,  The  Ohio  State  University,  Columbus, 
Ohio,  1988. 

3.  Hui-Huang  Abel  Chyou,  "Variational  Formulation  and  Finite  Element  Imple¬ 
mentation  of  Pagano's  Theory  of  Laminated  Plates”,  Ph.D.  Dissertation,  The 
Ohio  State  University,  Columbus,  Ohio,  1989. 

4.  Kazan  A.  Dandan,  ’’Finite  Element  Analysis  of  Laminated  Composite  Axisym- 
metric  Solids”,  M.S.  Thesis,  The  Ohio  State  University,  Columbus,  Ohio,  1988. 

5.  Mehdi  Moazzami,  ”A  Higher-Order  Layerwise  Theory  of  Laminated  Plates", 
Ph.D.  Dissertation,  The  Ohio  State  University,  Columbus,  Ohio,  1991. 


166 


☆  U  S.  GOVERNMENT  PRINTING  OFFICE:  IW  -  64S-U7/623M 


